# THE LAW BEHIND DISPUTE ONSET: 
# HOW LEGAL UNCERTAINTY DRIVES MARITIME BOUNDARY DISPUTES#

# Umut Yüksel

# Universitat Pompeu Fabra

# Contact: umut.yuksel@upf.edu or umut.yuksel@graduateinstitute.ch

# REPLICATION SCRIPT #
# Version 19 November 2024 # 

#### A. PACKAGES ####

library(alpaca)
library(brglm)
library(brglm2)
library(broom) 
library(clarify)
library(fixest)
library(foreign)
library(ggfixest)
library(ggpubr)
library(interactions)
library(jtools) 
library(knitr)
library(lme4)
library(lmtest)
library(logistf)
library(marginaleffects)
library(margins)
library(MASS)
library(mfx)
library(multiwayvcov)
library(plm)
library(prediction)
library(psych)
library(RColorBrewer)
library(readr) 
library(readxl)
library(sandwich) 
library(stargazer)
library(survival)
library(texreg)
library(tidyverse) 
library(vistime)
library(Zelig)

#### B. WORKING DIRECTORY ####

# Setting the working directory 
# (to wherever you have the following data files: "mbm_phi_jcr.RData" and "mbm_phi_wo_F0_jcr.RData")
setwd("~/TheLawBehind_Replication")

#### C. DATA ####

# Load full dataset:
load("mbm_phi_jcr.RData")
# Load dataset with fully delimited boundaries excluded:
load("mbm_phi_wo_F0_jcr.RData")

#### D. DESCRIPTIVES ####

# Number of unique dyads
mbm_phi |> group_by(dyad) |>
  count() 
# 444 x 2; 444 unique dyads

# Number of dyads per year:
no_dyads <- mbm_phi |> 
  group_by(year) |> 
  count(dyad) |> 
  summarise(no_dyads = sum(n))

# D.1. Statuses ####

# Counts of boundary statuses per year:
mbm_phi_statuses <- mbm_phi |> 
  group_by(year) |> 
  count(status_num)

# Merging number of dyads and statuses:
mbm_phi_statuses_and_numbers <- mbm_phi_statuses |> 
  right_join(no_dyads)

# Color palette:
my_greys_dist_plot = brewer.pal(n = 9, "Greys")[2:6]

# Status distribution:
status_distripution_plot <- ggplot(mbm_phi_statuses_and_numbers, 
                                   aes(x = as.numeric(year), 
                                       y = n, 
                                       fill = status_num)) +
  geom_bar(stat = "identity", position = "fill", width = 0.8) + 
  geom_line(aes(y = no_dyads/400), linetype = 1, linewidth = 0.75) + 
  theme_bw(base_size = 21) +
  xlab("Year") +
  scale_y_continuous(name = "Proportion of dyadic relations with a given status", 
                     breaks = seq(0, 1, 0.1),
                     sec.axis = sec_axis(transform = ~.*400,   
                                         name = "Number of dyads",
                                         breaks = seq(0, 400, 40))) + 
  theme(legend.position = "bottom",
        legend.justification = "center") +
  scale_x_continuous(breaks=seq(1946, 2016, 5)) +
  scale_fill_manual(values = my_greys_dist_plot) + 
  theme(legend.title = element_blank()) + 
  guides(fill = guide_legend(nrow = 3, 
                             byrow = TRUE, 
                             reverse = TRUE))

#### D.1.1. <<< Main FIGURE 1 >>> Status distribution ####
ggsave(status_distripution_plot, 
       filename = "main_figure_1.pdf", 
       width = 30, 
       height = 25, 
       units = "cm")

# Numbers:
status_numbers_8_years <- mbm_phi |>
  group_by(year, status_num) |>
  tally() |>
  filter(year %in% c(1946, 1956, 1966, 1976, 1986, 1996, 2006, 2016)) |>
  group_by(status_num) |>
  pivot_wider(names_from = year, values_from = n) 

# Proportions:
status_proportions_8_years <- mbm_phi |>
  group_by(year, status_num) |>
  tally() |>
  filter(year %in% c(1946, 1956, 1966, 1976, 1986, 1996, 2006, 2016)) |>
  group_by(year) |>
  add_tally(n) |>
  mutate(prop = round((n/nn), 2)) |>
  group_by(status_num) |>
  dplyr::select(year, status_num, prop) |>
  pivot_wider(names_from = year, values_from = prop) 

# Total number of dyads in dispute: 
dyads_in_dispute_numbers_8_years <- mbm_phi |>
  filter(year %in% c(1946, 1956, 1966, 1976, 1986, 1996, 2006, 2016)) |>
  filter(status_num %in% c("Undelimited & Disputed", "Partially Delimited & Disputed")) |>
  group_by(year) |>
  tally() 

# Total number of dyads: 
all_dyads_8_years <- mbm_phi |>
  filter(year %in% c(1946, 1956, 1966, 1976, 1986, 1996, 2006, 2016)) |>
  group_by(year) |>
  tally() 

#### D.1.2. <<< OS TABLE S3.1 >>> Status evolution, every 10 years ####
# The information contained in status_proportions_8_years,
# dyads_in_dispute_numbers_8_years, all_dyads_8_years are presented together in a table. 

# D.2. Number of onsets ####

# Number of onsets per year:
mbm_phi_ambilex_onsets <- mbm_phi |>
  group_by(year) |> 
  add_tally(onsmdisp, name = "number_of_onsets") |>
  group_by(year, number_of_onsets, ambilex) |>
  summarise(number_of_dyads = n_distinct(dyad)) |>
  mutate(prop_onsets = number_of_onsets/number_of_dyads)

# Save the annotations:
annotation <- data.frame(
  x = c((1946 + 1960)/2,
        (1960 + 1969)/2,
        (1969 + 1982)/2,
        (1982 + 1993)/2,
        (1993 + 2016)/2),
  y = rep(7.75, 5),
  label = c("LOW", "MEDIUM", "HIGH", "MEDIUM", "LOW"))

# Plotting onsets over time:
onset_plot <- ggplot(mbm_phi_ambilex_onsets, aes(x = year)) +
  geom_histogram(aes(y = number_of_onsets), 
                 stat = "identity", alpha = 0.75) +
  geom_line(aes(y = number_of_dyads/50), 
            linetype = 1, linewidth = 0.75) + 
  scale_y_continuous(name = "Number of disputes", 
                     breaks = seq(0, 8, 1), limits = c(0, 8),
                     sec.axis = sec_axis(trans = ~.*50,   
                                         name = "Number of dyads",
                                         breaks = seq(0, 400, 50))) + 
  scale_x_continuous(name = "Year of dispute onset", breaks = seq(1946, 2016, 5)) +
  theme_bw(base_size = 21) +
  geom_vline(xintercept = c(1960, 1993), linetype = 2, color = "red") + 
  geom_vline(xintercept = c(1969, 1982), linetype = 1, color = "red") +
  geom_label(data=annotation, aes(x = x, y = y, label = label), size = 6)

#### D.2.1. <<< Main FIGURE 2 >>> Number of onsets over time ####
# Save onsets over time:
ggsave(onset_plot, 
       filename = "main_figure_2.pdf",
       width = 30, 
       height = 15, units = "cm")

# D.3. Dispute onset descriptives ####

# How rare is this event?
mbm_phi |> summarize(mean(onsmdisp)) # 0.006105585 in the full dataset
mbm_phi_wo_F0 |> summarize(mean(onsmdisp)) # 0.007352941 if we exclude fully delimited boundries 

# How many dispute onsets? 
mbm_phi |> tally(onsmdisp) # 116

# How are disputes distributed across dyads? 
mbm_phi |> group_by(dyad) |>
  summarize(n_dispute = sum(onsmdisp, na.rm = TRUE)) |>
  filter(n_dispute > 0) |>
  arrange(desc(n_dispute)) |>
  group_by(n_dispute) |>
  count()
# 92 dyads with 1 dispute onsets
# 9 dyads with 2 dispute onsets
# 2 dyads with 3 dispute onsets

#### D.4. Summary statistics for variables used ####

#### D.4.1. <<< OS TABLE S3.2 >>> Descriptive Statistics ####
stargazer(mbm_phi_wo_F0, 
          keep = c("onsmdisp", 
                   "island_os","adjacency", "cseez",
                   "petrexpl", "petractive", 
                   "ntodel",
                   "amdisp", "atdisp", "joint_dem",
                   "cinc_ratio", "catch_tot_logged", "joint_unclos"),
          no.space = TRUE)
# This reports more than what is needed; the Table S3.2 focuses only on the listed variables. 

#### E. MODELS for BASELINE UNCERTAINTY ####

# Convert dyad to factor:
mbm_phi_wo_F0$dyad <- as.factor(mbm_phi_wo_F0$dyad)

#### E.1. FE Logit 1 ####

logit_fe1 <- alpaca::feglm(onsmdisp ~  ambilex_factor + 
                            ons_cntr + ons_cntr_sqd + ons_cntr_cbd | dyad | dyad, mbm_phi_wo_F0)

summary(logit_fe1)
logit_fe1_bc <- alpaca::biasCorr(logit_fe1) # Model
logit_fe1_bc_vcov <- vcov(logit_fe1_bc, type = "sandwich", cluster = ~ dyad)
logit_fe1_bc_vcov_rse <- sqrt(diag(logit_fe1_bc_vcov)) # Robust standard errors

summary_logit_fe1_bc_rcm <- summary(logit_fe1_bc, type = "sandwich", cluster = ~dyad)$cm
logit_fe1_bc_rpvalues <- summary_logit_fe1_bc_rcm[,4] # P-values

logit_fe1_bc_apes <- getAPEs(logit_fe1_bc)
logit_fe1_bc_apes_vcov <- vcov(logit_fe1_bc_apes, type = "sandwich", cluster = ~dyad)
logit_fe1_bc_apes_vcov_rse <- sqrt(diag(logit_fe1_bc_apes_vcov))

summary_logit_fe1_bc_apes_rcm <- summary(logit_fe1_bc_apes, type = "sandwich", cluster = ~dyad)$cm
logit_fe1_bc_apes_rpvalues <- summary_logit_fe1_bc_apes_rcm[,4]

# NB: Code for creating and presenting margins adapted from Supplementary material in 
# "Cook SJ, Hays JC, Franzese RJ. Fixed effects in rare events data: a penalized maximum likelihood solution. 
# Political Science Research and Methods. 2020;8(1):92-105. doi:10.1017/psrm.2018.40" 

mfx_logit_fe1 <- data.frame(Coefficient = summary_logit_fe1_bc_apes_rcm[,1],
                           SE = summary_logit_fe1_bc_apes_rcm[,2],
                           lower = summary_logit_fe1_bc_apes_rcm[,1] - summary_logit_fe1_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                           upper = summary_logit_fe1_bc_apes_rcm[,1] + summary_logit_fe1_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                           Model = "MFX Logit FE 1")
mfx_logit_fe1 <- data.frame(Variable_Contrast = row.names(mfx_logit_fe1), mfx_logit_fe1)
rownames(mfx_logit_fe1) = NULL

#### E.2. FE Logit 2 ####

logit_fe2 <- alpaca::feglm(onsmdisp ~ ambilex_factor + 
                            ntodel_lag_log + 
                            unclos_status + 
                            petrexpl + 
                            petractive_lag +
                            amdisp_lag + 
                            atdisp_lag +
                            joint_dem +
                            ons_cntr + ons_cntr_sqd + ons_cntr_cbd | dyad | dyad, mbm_phi_wo_F0)
summary(logit_fe2)
logit_fe2_bc <- alpaca::biasCorr(logit_fe2)

logit_fe2_bc_vcov <- vcov(logit_fe2_bc, type = "sandwich", cluster = ~dyad)
logit_fe2_bc_vcov_rse <- sqrt(diag(logit_fe2_bc_vcov))

summary_logit_fe2_bc_rcm <- summary(logit_fe2_bc, type = "sandwich", cluster = ~dyad)$cm
logit_fe2_bc_rpvalues <- summary_logit_fe2_bc_rcm[,4]

logit_fe2_bc_apes <- getAPEs(logit_fe2_bc)

logit_fe2_bc_apes_vcov <- vcov(logit_fe2_bc_apes, type = "sandwich", cluster = ~dyad)
logit_fe2_bc_apes_vcov_rse <- sqrt(diag(logit_fe2_bc_apes_vcov))

summary_logit_fe2_bc_apes_rcm <- summary(logit_fe2_bc_apes, type = "sandwich", cluster = ~dyad)$cm
logit_fe2_bc_apes_rpvalues <- summary_logit_fe2_bc_apes_rcm[,4]

mfx_logit_fe2 <- data.frame(Coefficient = summary_logit_fe2_bc_apes_rcm[,1],
                           SE = summary_logit_fe2_bc_apes_rcm[,2],
                           lower = summary_logit_fe2_bc_apes_rcm[,1] - summary_logit_fe2_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                           upper = summary_logit_fe2_bc_apes_rcm[,1] + summary_logit_fe2_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                           Model = "MFX Logit FE 2")
mfx_logit_fe2 <- data.frame(Variable_Contrast = row.names(mfx_logit_fe2), mfx_logit_fe2)
rownames(mfx_logit_fe2) = NULL

Variable_Contrast_logit_fe <- mfx_logit_fe2$Variable_Contrast

#### E.3. FE Rare Logit 1 ####

# Model: 

rare_logit_fe1 <- glm(onsmdisp ~ ambilex_factor +
                        ons_cntr +
                        ons_cntr_sqd +
                        ons_cntr_cbd + 
                        factor(dyad),
                      data = mbm_phi_wo_F0, 
                      family = binomial(link = "logit"), 
                      method = "brglmFit")

# Get summary:
coeftest(rare_logit_fe1, vcov. = vcovCL(rare_logit_fe1, cluster = mbm_phi_wo_F0$dyad))

# Save robust standard errors:
cov_rare_logit_fe1 <- vcovCL(rare_logit_fe1, cluster = mbm_phi_wo_F0$dyad)

robust_se_rare_logit_fe1 <- sqrt(diag(cov_rare_logit_fe1))

# Margins:
mfx_rare_logit_fe1 <- margins(rare_logit_fe1, 
                              variables = c("ambilex_factor",
                                            "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                              vcov = cov_rare_logit_fe1)

mfx_rare_logit_fe1_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_fe1)[,1],
                                    Coefficient = summary(mfx_rare_logit_fe1)[,2],
                                    SE = summary(mfx_rare_logit_fe1)[,3],
                                    lower = summary(mfx_rare_logit_fe1)[,6],
                                    upper = summary(mfx_rare_logit_fe1)[,7],
                                    Model = "MFX Rare Logit FE 1")

rownames(mfx_rare_logit_fe1_df) <- NULL

#### E.4. FE Rare Logit 2 ####

# Model: 
rare_logit_fe2 <- glm(onsmdisp ~ ambilex_factor + 
                        ntodel_lag_log + 
                        unclos_status + 
                        petrexpl + 
                        petractive_lag +
                        amdisp_lag + 
                        atdisp_lag + 
                        joint_dem +
                        ons_cntr + 
                        ons_cntr_sqd + 
                        ons_cntr_cbd + 
                        factor(dyad),
                      data = mbm_phi_wo_F0, 
                      family = binomial(link = "logit"), 
                      method = "brglmFit")

# Get summary:
coeftest(rare_logit_fe2, vcov. = vcovCL(rare_logit_fe2, 
                                        cluster = mbm_phi_wo_F0$dyad))

# Save robust standard errors:
cov_rare_logit_fe2 <- vcovCL(rare_logit_fe2, 
                             cluster = mbm_phi_wo_F0$dyad)

robust_se_rare_logit_fe2 <- sqrt(diag(cov_rare_logit_fe2))

# Margins:
mfx_rare_logit_fe2 <- margins(rare_logit_fe2, 
                              variables = c("ambilex_factor",
                                            "ntodel_lag_log",
                                            "unclos_status", 
                                            "petrexpl", "petractive_lag",
                                            "amdisp_lag", "atdisp_lag", "joint_dem",
                                            "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                              vcov = cov_rare_logit_fe2)

mfx_rare_logit_fe2_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_fe2)[,1],
                                    Coefficient = summary(mfx_rare_logit_fe2)[,2],
                                    SE = summary(mfx_rare_logit_fe2)[,3],
                                    lower = summary(mfx_rare_logit_fe2)[,6],
                                    upper = summary(mfx_rare_logit_fe2)[,7],
                                    Model = "MFX Rare Logit FE 2")

rownames(mfx_rare_logit_fe2_df) <- NULL

#### E.5. Regresion output for models E.1 - E.4 #### 

# Selecting GOFs for the rare event logit models:
rare_logit_fe1_sle <- texreg::extract(rare_logit_fe1)
rare_logit_fe1_sle@gof.names <- rare_logit_fe1_sle@gof.names[-(1:8)]
rare_logit_fe1_sle@gof.decimal <- rare_logit_fe1_sle@gof.decimal[-(1:8)]
rare_logit_fe1_sle@gof <- rare_logit_fe1_sle@gof[-(1:8)]

rare_logit_fe2_sle <- texreg::extract(rare_logit_fe2)
rare_logit_fe2_sle@gof.names <- rare_logit_fe2_sle@gof.names[-(1:8)]
rare_logit_fe2_sle@gof.decimal <- rare_logit_fe2_sle@gof.decimal[-(1:8)]
rare_logit_fe2_sle@gof <- rare_logit_fe2_sle@gof[-(1:8)]

#### E.5.1. <<< Main TABLE 2 >>> Regression output for models E.1 - E.4 #### 

texreg(list(logit_fe1_bc, rare_logit_fe1_sle, logit_fe2_bc, rare_logit_fe2_sle), 
       override.se = list(logit_fe1_bc_vcov_rse, robust_se_rare_logit_fe1, logit_fe2_bc_vcov_rse, robust_se_rare_logit_fe2), 
       override.pvalues = list(logit_fe1_bc_rpvalues, NULL, logit_fe2_bc_rpvalues, NULL),
       # single.row = TRUE,
       no.margin = TRUE,
       fontsize = "scriptsize",
       custom.model.names = c("FE Logit", "Penalized ML", "FE Logit", "Penalized ML"),
       custom.coef.map = list("ambilex_factorMedium" = "MEDIUM legal uncertainty (ref: LOW)",
                              "ambilex_factorHigh" = "HIGH legal uncertainty (ref: LOW)",
                              "ntodel_lag_log" = "No. of dyadic boundaries to delimit (lagged & logged)",
                              "unclos_statusJoint UNCLOS" = "Joint UNCLOS (ref: Before UNCLOS)", 
                              "unclos_statusNeither UNCLOS" = "Neither UNCLOS (ref: Before UNCLOS)", 
                              "unclos_statusOne state UNCLOS" = "One state UNCLOS (ref: Before UNCLOS)", 
                              "petrexpl" = "Hydrocarbon exploration",
                              "petractive_lag" = "Hydrocarbon activity",
                              "amdisp_lag" = "Ongoing maritime boundary dispute (lagged)", 
                              "atdisp_lag" = "Ongoing related territorial dispute (lagged)",
                              "joint_demJoint democracy" = "One side DEMOCRATIC (ref: joint AUTOCRACY)", 
                              "joint_demOne side democratic" = "Joint DEMOCRACY (ref: joint AUTOCRACY)"), 
       # omit.coef = "(dyad)|(ons_cntr)|(ons_cntr_cbd)|(ons_cntr_sqd)",
       # reorder.coef = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13),
       # gof.names = list("", "", "", "", "", "", ""),
       custom.gof.rows = list("Dyad-fixed effects" = c("Yes", "Yes", "Yes", "Yes"),
                              "Year-fixed effects" = c("No", "No", "No", "No"),
                              "Cubic polynomial for time" = c("Yes", "Yes", "Yes", "Yes")),
       include.deviance = FALSE,
       booktabs = TRUE,
       dcolumn = TRUE,
       stars = c(0.001, 0.01, 0.05, 0.10),
       symbol = "+")

#### E.6. Marginal effect plot for models E.1 - E.4 #### 

# NB: Code for creating and presenting margins adapted from Supplementary material in 
# "Cook SJ, Hays JC, Franzese RJ. Fixed effects in rare events data: a penalized maximum likelihood solution. 
# Political Science Research and Methods. 2020;8(1):92-105. doi:10.1017/psrm.2018.40" 

# Margins logit_fe and rare events logits:
margins_logit_fe_rare_logit_fe <- data.frame(rbind(mfx_logit_fe1, mfx_rare_logit_fe1_df,
                                              mfx_logit_fe2, mfx_rare_logit_fe2_df)) |>
  filter(Variable_Contrast != "ons_cntr" & 
           Variable_Contrast != "ons_cntr_sqd" & 
           Variable_Contrast != "ons_cntr_cbd")

interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

margins_logit_fe_rare_logit_fe$Variable_Contrast <- factor(margins_logit_fe_rare_logit_fe$Variable_Contrast, 
                                                      levels = c("ambilex_factorMedium", 
                                                                 "ambilex_factorHigh", 
                                                                 "ntodel_lag_log", 
                                                                 "unclos_statusJoint UNCLOS", 
                                                                 "unclos_statusNeither UNCLOS", 
                                                                 "unclos_statusOne state UNCLOS",
                                                                 "petrexpl", 
                                                                 "petractive_lag", 
                                                                 "amdisp_lag", 
                                                                 "atdisp_lag", 
                                                                 "joint_demOne side democratic",
                                                                 "joint_demJoint democracy"),
                                                      labels = c("MEDIUM legal uncertainty (ref: LOW)",
                                                                 "HIGH legal uncertainty (ref: LOW)",
                                                                 "No. of dyadic boundaries to delimit (lagged & logged)",
                                                                 "Joint UNCLOS (ref: Before UNCLOS)", 
                                                                 "Neither UNCLOS (ref: Before UNCLOS)", 
                                                                 "One state UNCLOS (ref: Before UNCLOS)", 
                                                                 "Hydrocarbon exploration",
                                                                 "Hydrocarbon activity",
                                                                 "Ongoing maritime boundary dispute (lagged)", 
                                                                 "Ongoing related territorial dispute (lagged)",
                                                                 "One side DEMOCRATIC (ref: joint AUTOCRACY)", 
                                                                 "Joint DEMOCRACY (ref: joint AUTOCRACY)"))

margins_logit_fe_rare_logit_fe$Model <- factor(margins_logit_fe_rare_logit_fe$Model, 
                                          levels = c("MFX Rare Logit FE 2", 
                                                     "MFX Rare Logit FE 1", 
                                                     "MFX Logit FE 2",
                                                     "MFX Logit FE 1"), 
                                          labels = c("(4) Penalized ML FE with controls",
                                                     "(3) Penalized ML FE" , 
                                                     "(2) Logit FE with controls",
                                                     "(1) Logit FE"))
# Defining tones of grey: 
my_greys = brewer.pal(n = 9, "Greys")[6:9] # excluding two lighter hues

mep_main_all <- ggplot(margins_logit_fe_rare_logit_fe, aes(colour = Model)) + 
  scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1,
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Model),
                  lwd = 1/2, position = position_dodge(width = 0.5), 
                  fill = "white", size = 1) + 
  scale_shape(solid = TRUE) + scale_x_discrete(limits = rev(levels(margins_logit_fe_rare_logit_fe$Variable_Contrast))) + 
  theme_bw(base_size = 30) +
  coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 2)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 2)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_main_all

#### E.6.1. <<< OS FIGURE S4.1 >>> All margins for models with baseline uncertainty ####
# Saving: 
ggsave(mep_main_all, 
       filename = "os_figure_S4_1.pdf",
       width = 50, 
       height = 40, units = "cm")

# Only IVs: 

# Margins logit and rare events logits:
margins_logit_fe_rare_logit_fe <- data.frame(rbind(mfx_logit_fe1, mfx_rare_logit_fe1_df,
                                              mfx_logit_fe2, mfx_rare_logit_fe2_df)) |>
  filter(Variable_Contrast == "ambilex_factorMedium" | 
           Variable_Contrast == "ambilex_factorHigh")

margins_logit_fe_rare_logit_fe$Variable_Contrast <- factor(margins_logit_fe_rare_logit_fe$Variable_Contrast, 
                                                      levels = c("ambilex_factorMedium", 
                                                                 "ambilex_factorHigh"),
                                                      labels = c("MEDIUM legal uncertainty (ref: LOW)",
                                                                 "HIGH legal uncertainty (ref: LOW)"))

margins_logit_fe_rare_logit_fe$Model <- factor(margins_logit_fe_rare_logit_fe$Model, 
                                          levels = c("MFX Rare Logit FE 2", 
                                                     "MFX Rare Logit FE 1", 
                                                     "MFX Logit FE 2",
                                                     "MFX Logit FE 1"), 
                                          labels = c("(4) Penalized ML FE with controls",
                                                     "(3) Penalized ML FE" , 
                                                     "(2) Logit FE with controls",
                                                     "(1) Logit FE"))

mep_main_ambilex <- ggplot(margins_logit_fe_rare_logit_fe, aes(colour = Model)) + scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1,
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, 
                 position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Model), 
                  lwd = 1/2, 
                  position = position_dodge(width = 0.5), 
                  fill = "white", 
                  size = 1) + 
  scale_shape(solid = TRUE)+ 
  scale_x_discrete(limits = rev(levels(margins_logit_fe_rare_logit_fe$Variable_Contrast))) + 
  theme_bw(base_size = 30) + 
  coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 2)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 2)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_main_ambilex

#### E.6.2. <<< Main FIGURE 3 >>> Uncertainty related variables ####
# Saving: 
ggsave(mep_main_ambilex, 
       filename = "main_figure_3.pdf",
       width = 50, 
       height = 20, units = "cm")

#### F. MODELS for DYADIC UNCERTAINTY  ####

#### F.1. Year FE Logit 1 ####

year_fes1 <- alpaca::feglm(onsmdisp ~ island_os + 
                             adjacency + 
                             cseez + 
                             ons_cntr +
                             ons_cntr_sqd + 
                             ons_cntr_cbd | year | year, mbm_phi_wo_F0) 

summary(year_fes1)

year_fes1_bc <- alpaca::biasCorr(year_fes1) # Model
year_fes1_bc_vcov <- vcov(year_fes1_bc, type = "sandwich", cluster = ~ year)
year_fes1_bc_vcov_rse <- sqrt(diag(year_fes1_bc_vcov)) # Robust standard errors

summary_year_fes1_bc_rcm <- summary(year_fes1_bc, type = "sandwich", cluster = ~ year)$cm
year_fes1_bc_rpvalues <- summary_year_fes1_bc_rcm[,4] # P-values

year_fes1_bc_apes <- getAPEs(year_fes1_bc)
year_fes1_bc_apes_vcov <- vcov(year_fes1_bc_apes, type = "sandwich", cluster = ~ year)
year_fes1_bc_apes_vcov_rse <- sqrt(diag(year_fes1_bc_apes_vcov))

summary_year_fes1_bc_apes_rcm <- summary(year_fes1_bc_apes, type = "sandwich", cluster = ~ year)$cm
year_fes1_bc_apes_rpvalues <- summary_year_fes1_bc_apes_rcm[,4]

mfx_year_fes1 <- data.frame(Coefficient = summary_year_fes1_bc_apes_rcm[,1],
                            SE = summary_year_fes1_bc_apes_rcm[,2],
                            lower = summary_year_fes1_bc_apes_rcm[,1] - summary_year_fes1_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                            upper = summary_year_fes1_bc_apes_rcm[,1] + summary_year_fes1_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                            Model = "MFX Year FE 1")

mfx_year_fes1 <- data.frame(Variable_Contrast = row.names(mfx_year_fes1), mfx_year_fes1)
rownames(mfx_year_fes1) = NULL

#### F.2. Year FE Logit 2 ####

year_fes2 <- alpaca::feglm(onsmdisp ~ island_os + 
                             adjacency + 
                             cseez + 
                             ntodel_lag_log + 
                             petrexpl + 
                             petractive_lag +
                             catch_tot_lag_logged + 
                             amdisp_lag + 
                             atdisp_lag +
                             joint_dem +
                             cinc_ratio +
                             ons_cntr +
                             ons_cntr_sqd + 
                             ons_cntr_cbd | year | year, mbm_phi_wo_F0) 

summary(year_fes2)

year_fes2_bc <- alpaca::biasCorr(year_fes2) # Model
year_fes2_bc_vcov <- vcov(year_fes2_bc, type = "sandwich", cluster = ~ year)
year_fes2_bc_vcov_rse <- sqrt(diag(year_fes2_bc_vcov)) # Robust standard errors

summary_year_fes2_bc_rcm <- summary(year_fes2_bc, type = "sandwich", cluster = ~ year)$cm
year_fes2_bc_rpvalues <- summary_year_fes2_bc_rcm[,4] # P-values

year_fes2_bc_apes <- getAPEs(year_fes2_bc)
year_fes2_bc_apes_vcov <- vcov(year_fes2_bc_apes, type = "sandwich", cluster = ~ year)
year_fes2_bc_apes_vcov_rse <- sqrt(diag(year_fes2_bc_apes_vcov))

summary_year_fes2_bc_apes_rcm <- summary(year_fes2_bc_apes, type = "sandwich", cluster = ~ year)$cm
year_fes2_bc_apes_rpvalues <- summary_year_fes2_bc_apes_rcm[,4]

mfx_year_fes2 <- data.frame(Coefficient = summary_year_fes2_bc_apes_rcm[,1],
                            SE = summary_year_fes2_bc_apes_rcm[,2],
                            lower = summary_year_fes2_bc_apes_rcm[,1] - summary_year_fes2_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                            upper = summary_year_fes2_bc_apes_rcm[,1] + summary_year_fes2_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                            Model = "MFX Year FE 2")

mfx_year_fes2 <- data.frame(Variable_Contrast = row.names(mfx_year_fes2), mfx_year_fes2)
rownames(mfx_year_fes2) = NULL

#### F.3. Year FE Rare Logit 1 ####

# Model: 
rare_logit_year_fes1 <- glm(onsmdisp ~ island_os + 
                              adjacency + 
                              cseez + 
                              ons_cntr +
                              ons_cntr_sqd +
                              ons_cntr_cbd + 
                              factor(year),
                            data = mbm_phi_wo_F0, 
                            family = binomial(link = "logit"), 
                            method = "brglmFit")

# Get summary:
coeftest(rare_logit_year_fes1, vcov. = vcovCL(rare_logit_year_fes1, cluster = mbm_phi_wo_F0$year))

# Save robust standard errors:
cov_rare_logit_year_fes1 <- vcovCL(rare_logit_year_fes1, 
                                   cluster = mbm_phi_wo_F0$year)

robust_se_rare_logit_year_fes1 <- sqrt(diag(cov_rare_logit_year_fes1))

# Margins:
mfx_rare_logit_year_fes1 <- margins(rare_logit_year_fes1, 
                                    variables = c("island_os", "adjacency", "cseez", 
                                                  "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                                    vcov = cov_rare_logit_year_fes1)

mfx_rare_logit_year_fes1_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_year_fes1)[,1],
                                          Coefficient = summary(mfx_rare_logit_year_fes1)[,2],
                                          SE = summary(mfx_rare_logit_year_fes1)[,3],
                                          lower = summary(mfx_rare_logit_year_fes1)[,6],
                                          upper = summary(mfx_rare_logit_year_fes1)[,7],
                                          Model = "MFX Rare Logit Year FE 1")

rownames(mfx_rare_logit_year_fes1_df) <- NULL

#### F.4. Year FE Rare Logit 2 ####

# Model: 
rare_logit_year_fes2 <- glm(onsmdisp ~ island_os + 
                              adjacency + 
                              cseez + 
                              ntodel_lag_log + 
                              # unclos_status + 
                              petrexpl + 
                              petractive_lag +
                              catch_tot_lag_logged + 
                              amdisp_lag + 
                              atdisp_lag +
                              joint_dem +
                              cinc_ratio +
                              ons_cntr +
                              ons_cntr_sqd + 
                              ons_cntr_cbd +
                              factor(year),
                            data = mbm_phi_wo_F0, 
                            family = binomial(link = "logit"), 
                            method = "brglmFit", 
                            pl = TRUE)

# Get summary:
coeftest(rare_logit_year_fes2, vcov. = vcovCL(rare_logit_year_fes2, cluster = mbm_phi_wo_F0$year))

# Save robust standard errors:
cov_rare_logit_year_fes2 <- vcovCL(rare_logit_year_fes2, 
                                   cluster = mbm_phi_wo_F0$year)

robust_se_rare_logit_year_fes2 <- sqrt(diag(cov_rare_logit_year_fes2))

# Margins:
mfx_rare_logit_year_fes2 <- margins(rare_logit_year_fes2, 
                                    variables = c("island_os", "adjacency", "cseez", 
                                                  "ntodel_lag_log", "petrexpl", "petractive_lag",
                                                  "catch_tot_lag_logged", "amdisp_lag", "atdisp_lag",
                                                  "joint_dem", "cinc_ratio",
                                                  "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                                    vcov = cov_rare_logit_year_fes2)

mfx_rare_logit_year_fes2_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_year_fes2)[,1],
                                          Coefficient = summary(mfx_rare_logit_year_fes2)[,2],
                                          SE = summary(mfx_rare_logit_year_fes2)[,3],
                                          lower = summary(mfx_rare_logit_year_fes2)[,6],
                                          upper = summary(mfx_rare_logit_year_fes2)[,7],
                                          Model = "MFX Rare Logit Year FE 2")

rownames(mfx_rare_logit_year_fes2_df) <- NULL

#### F.5. Regression output for models E.1 - E.4 ####

# Selecting GOFs for the rare event logit models:
rare_logit_year_fes1_sle <- texreg::extract(rare_logit_year_fes1)
rare_logit_year_fes1_sle@gof.names <- rare_logit_year_fes1_sle@gof.names[-(1:8)]
rare_logit_year_fes1_sle@gof.decimal <- rare_logit_year_fes1_sle@gof.decimal[-(1:8)]
rare_logit_year_fes1_sle@gof <- rare_logit_year_fes1_sle@gof[-(1:8)]

rare_logit_year_fes2_sle <- texreg::extract(rare_logit_year_fes2)
rare_logit_year_fes2_sle@gof.names <- rare_logit_year_fes2_sle@gof.names[-(1:8)]
rare_logit_year_fes2_sle@gof.decimal <- rare_logit_year_fes2_sle@gof.decimal[-(1:8)]
rare_logit_year_fes2_sle@gof <- rare_logit_year_fes2_sle@gof[-(1:8)]

#### F.5.1 <<< Main TABLE 3 >>> Regression output for models E.1 - E.4 ####

texreg(list(year_fes1_bc, rare_logit_year_fes1_sle, year_fes2_bc, rare_logit_year_fes2_sle),
       override.se = list(year_fes1_bc_vcov_rse, robust_se_rare_logit_year_fes1, year_fes2_bc_vcov_rse, robust_se_rare_logit_year_fes2), 
       override.pvalues = list(year_fes1_bc_rpvalues, NULL, year_fes2_bc_rpvalues, NULL),
       # single.row = TRUE,
       no.margin = TRUE,
       fontsize = "scriptsize",
       custom.model.names = c("FE Logit", "Penalized ML", "FE Logit", "Penalized ML"),
       custom.coef.map = list("island_os" = "Offshore island",
                              "adjacency" = "Adjacency of coasts",
                              "cseez" = "Beyond territorial sea", 
                              "ntodel_lag_log" = "No. of dyadic boundaries to delimit (lagged & logged)",
                              "petrexpl" = "Hydrocarbon exploration",
                              "petractive_lag" = "Hydrocarbon activity",
                              "catch_tot_lag_logged" = "Total catch (lagged & logged)",
                              "cinc_ratio" = "Capability ratio",
                              "amdisp_lag" = "Ongoing maritime boundary dispute (lagged)", 
                              "atdisp_lag" = "Ongoing related territorial dispute (lagged)",
                              "joint_demJoint democracy" = "One side DEMOCRATIC (ref: joint AUTOCRACY)", 
                              "joint_demOne side democratic" = "Joint DEMOCRACY (ref: joint AUTOCRACY)"), 
       # omit.coef = "(dyad)|(ons_cntr)|(ons_cntr_cbd)|(ons_cntr_sqd)",
       # reorder.coef = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13),
       gof.names = list("", "", "", "", "", "", ""),
       custom.gof.rows = list("Dyad-fixed effects" = c("No", "No", "No", "No"),
                              "Year-fixed effects" = c("Yes", "Yes", "Yes", "Yes"),
                              "Cubic polynomial for time" = c("Yes", "Yes", "Yes", "Yes")),
       include.deviance = FALSE,
       booktabs = TRUE,
       dcolumn = TRUE, 
       stars = c(0.001, 0.01, 0.05, 0.10),
       symbol = "+")

#### F.6. Marginal effect plot for models E.1 - E.4 ####

# Margins fe logit and rare events logit:
margins_logit_rare_logit_year_fe <- data.frame(rbind(mfx_year_fes1, mfx_rare_logit_year_fes1_df,
                                                     mfx_year_fes2, mfx_rare_logit_year_fes2_df)) |>
  filter(Variable_Contrast != "ons_cntr" & 
           Variable_Contrast != "ons_cntr_sqd" & 
           Variable_Contrast != "ons_cntr_cbd")

interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

margins_logit_rare_logit_year_fe$Variable_Contrast <- factor(margins_logit_rare_logit_year_fe$Variable_Contrast, 
                                                           levels = c("island_os",
                                                                      "adjacency",
                                                                      "cseez",
                                                                      "ntodel_lag_log", 
                                                                      "petrexpl", 
                                                                      "petractive_lag", 
                                                                      "amdisp_lag", 
                                                                      "atdisp_lag", 
                                                                      "joint_demOne side democratic",
                                                                      "joint_demJoint democracy",
                                                                      "catch_tot_lag_logged", 
                                                                      "cinc_ratio"),
                                                           labels = c("Dyadic uncert.: Offshore island",
                                                                      "Dyadic uncert.: Adjacency of coasts", 
                                                                      "Dyadic uncert.: Beyond territorial sea", 
                                                                      "No. of dyadic boundaries to delimit (lagged & logged)",
                                                                      "Hydrocarbon exploration",
                                                                      "Hydrocarbon activity",
                                                                      "Ongoing maritime boundary dispute (lagged)", 
                                                                      "Ongoing related territorial dispute (lagged)",
                                                                      "One side DEMOCRATIC (ref: joint AUTOCRACY)", 
                                                                      "Joint DEMOCRACY (ref: joint AUTOCRACY)",
                                                                      "Total catch (lagged & logged)",
                                                                      "Capability ratio"))

margins_logit_rare_logit_year_fe$Model <- factor(margins_logit_rare_logit_year_fe$Model, 
                                               levels = c("MFX Rare Logit Year FE 2", 
                                                          "MFX Rare Logit Year FE 1", 
                                                          "MFX Year FE 2",
                                                          "MFX Year FE 1"), 
                                               labels = c("(4) Penalized ML FE with controls",
                                                          "(3) Penalized ML FE" , 
                                                          "(2) Logit FE with controls",
                                                          "(1) Logit FE"))

mep_dyadic_all <- ggplot(margins_logit_rare_logit_year_fe, aes(colour = Model)) + 
  scale_colour_manual(values = my_greys) +
  geom_hline(yintercept = 0, colour = "red", lty = 2) + geom_linerange(aes(x = Variable_Contrast,
                                ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                                lwd = 1,
                                position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Model), 
                  lwd = 1/2, position = position_dodge(width = 0.5), 
                  fill = "white", size = 1) + 
  scale_shape(solid = TRUE) + scale_x_discrete(limits = rev(levels(margins_logit_rare_logit_year_fe$Variable_Contrast))) + 
  theme_bw(base_size = 30) + coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 2)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 2)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_dyadic_all

#### F.6.1. <<< OS FIGURE F.4.2 >>> All margins, models with dyadic uncertainty ####
ggsave(mep_dyadic_all, 
       filename = "os_figure_S4_2.pdf",
       width = 50, 
       height = 40, units = "cm")

# Only ambilex:

# Margins logit fe and rare events logits:
margins_logit_rare_logit_year_fe <- data.frame(rbind(mfx_year_fes1, mfx_rare_logit_year_fes1_df,
                                                   mfx_year_fes2, mfx_rare_logit_year_fes2_df)) |>
  filter(Variable_Contrast == "island_os" | Variable_Contrast == "adjacency" |
           Variable_Contrast == "cseez")

interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

margins_logit_rare_logit_year_fe$Variable_Contrast <- factor(margins_logit_rare_logit_year_fe$Variable_Contrast, 
                                                           levels = c("island_os",
                                                                      "adjacency",
                                                                      "cseez"),
                                                           labels = c("Dyadic uncert.: Offshore island",
                                                                      "Dyadic uncert.: Adjacency of coasts", 
                                                                      "Dyadic uncert.: Beyond territorial sea"))

margins_logit_rare_logit_year_fe$Model <- factor(margins_logit_rare_logit_year_fe$Model, 
                                               levels = c("MFX Rare Logit Year FE 2", 
                                                          "MFX Rare Logit Year FE 1", 
                                                          "MFX Year FE 2",
                                                          "MFX Year FE 1"), 
                                               labels = c("(4) Penalized ML FE with controls",
                                                          "(3) Penalized ML FE" , 
                                                          "(2) Logit FE with controls",
                                                          "(1) Logit FE"))

mep_dyadic_ambilex <- ggplot(margins_logit_rare_logit_year_fe, aes(colour = Model)) + scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1,
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Model), 
                  lwd = 1/2, position = position_dodge(width = 0.5), 
                  fill = "white", size = 1) + 
  scale_shape(solid = TRUE) + scale_x_discrete(limits = rev(levels(margins_logit_rare_logit_year_fe$Variable_Contrast))) + 
  theme_bw(base_size = 30) +  
  coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 2)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 2)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_dyadic_ambilex

#### F.6.2. <<< Main FIGURE 4 >>> Uncertainty related variables ####
ggsave(mep_dyadic_ambilex, 
       filename = "main_figure_4.pdf",
       width = 50, 
       height = 25, units = "cm")

#### G. ROBUSTNESS CHECKS ####

#### G.1. Other period variables ####

mbm_phi_wo_F0$same_legal_system[is.na(mbm_phi_wo_F0$same_legal_system)] = 0 # change NA to 0

mbm_phi_wo_F0 <- mbm_phi_wo_F0 %>% 
  mutate(unclos_signed = if_else(year >= 1982, 1, 0)) |>
  mutate(unclos_in_force = if_else(year >= 1995, 1, 0)) |> 
  mutate(joint_unclos = if_else(unclos_status == "Joint UNCLOS", 1, 0)) |>
  mutate(post_cw = if_else(year > 1991, 1, 0)) 

# G.1.1. UNCLOS negotiations + dyad FE ####

logit_period1 <- alpaca::feglm(onsmdisp ~ ambilex_factor + 
                                 ntodel_lag_log + 
                                 unclos_status + 
                                 petrexpl + 
                                 petractive_lag +
                                 amdisp_lag + 
                                 atdisp_lag +
                                 joint_dem +
                                 unclos_all_neg + 
                                 ons_cntr + ons_cntr_sqd + ons_cntr_cbd | dyad | dyad, mbm_phi_wo_F0)
summary(logit_period1)
logit_period1_bc <- alpaca::biasCorr(logit_period1)

logit_period1_bc_vcov <- vcov(logit_period1_bc, type = "sandwich", cluster = ~dyad)
logit_period1_bc_vcov_rse <- sqrt(diag(logit_period1_bc_vcov))

summary_logit_period1_bc_rcm <- summary(logit_period1_bc, type = "sandwich", cluster = ~dyad)$cm
logit_period1_bc_rpvalues <- summary_logit_period1_bc_rcm[,4]

logit_period1_bc_apes <- getAPEs(logit_period1_bc)

logit_period1_bc_apes_vcov <- vcov(logit_period1_bc_apes, type = "sandwich", cluster = ~dyad)
logit_period1_bc_apes_vcov_rse <- sqrt(diag(logit_period1_bc_apes_vcov))

summary_logit_period1_bc_apes_rcm <- summary(logit_period1_bc_apes, type = "sandwich", cluster = ~dyad)$cm
logit_period1_bc_apes_rpvalues <- summary_logit_period1_bc_apes_rcm[,4]

mfx_logit_period1 <- data.frame(Coefficient = summary_logit_period1_bc_apes_rcm[,1],
                              SE = summary_logit_period1_bc_apes_rcm[,2],
                              lower = summary_logit_period1_bc_apes_rcm[,1] - summary_logit_period1_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                              upper = summary_logit_period1_bc_apes_rcm[,1] + summary_logit_period1_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                              Model = "MFX Period Logit 1")
mfx_logit_period1 <- data.frame(Variable_Contrast = row.names(mfx_logit_period1), mfx_logit_period1)
rownames(mfx_logit_period1) = NULL

Variable_Contrast_logit <- mfx_logit_period1$Variable_Contrast

# G.1.2. UNCLOS in force + dyad FE ####

logit_period2 <- alpaca::feglm(onsmdisp ~ ambilex_factor +
                                 ntodel_lag_log + 
                                 petrexpl + 
                                 petractive_lag +
                                 amdisp_lag + 
                                 atdisp_lag +
                                 joint_dem +
                                 unclos_signed + 
                                 unclos_in_force +
                                 ons_cntr + ons_cntr_sqd + ons_cntr_cbd | dyad | dyad, mbm_phi_wo_F0)
summary(logit_period2)
logit_period2_bc <- alpaca::biasCorr(logit_period2)

logit_period2_bc_vcov <- vcov(logit_period2_bc, type = "sandwich", cluster = ~dyad)
logit_period2_bc_vcov_rse <- sqrt(diag(logit_period2_bc_vcov))

summary_logit_period2_bc_rcm <- summary(logit_period2_bc, type = "sandwich", cluster = ~dyad)$cm
logit_period2_bc_rpvalues <- summary_logit_period2_bc_rcm[,4]

logit_period2_bc_apes <- getAPEs(logit_period2_bc)

logit_period2_bc_apes_vcov <- vcov(logit_period2_bc_apes, type = "sandwich", cluster = ~dyad)
logit_period2_bc_apes_vcov_rse <- sqrt(diag(logit_period2_bc_apes_vcov))

summary_logit_period2_bc_apes_rcm <- summary(logit_period2_bc_apes, type = "sandwich", cluster = ~dyad)$cm
logit_period2_bc_apes_rpvalues <- summary_logit_period2_bc_apes_rcm[,4]

mfx_logit_period2 <- data.frame(Coefficient = summary_logit_period2_bc_apes_rcm[,1],
                              SE = summary_logit_period2_bc_apes_rcm[,2],
                              lower = summary_logit_period2_bc_apes_rcm[,1] - summary_logit_period2_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                              upper = summary_logit_period2_bc_apes_rcm[,1] + summary_logit_period2_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                              Model = "MFX Period Logit 2")
mfx_logit_period2 <- data.frame(Variable_Contrast = row.names(mfx_logit_period2), mfx_logit_period2)
rownames(mfx_logit_period2) = NULL

#### G.1.3. PostCW + dyad FE ####

logit_period3 <- alpaca::feglm(onsmdisp ~ ambilex_factor + 
                               ntodel_lag_log + 
                               unclos_status + 
                               petrexpl + 
                               petractive_lag +
                               amdisp_lag + 
                               atdisp_lag +
                               joint_dem +
                               post_cw +
                               ons_cntr + 
                               ons_cntr_sqd + 
                               ons_cntr_cbd | dyad | dyad, mbm_phi_wo_F0)
summary(logit_period3)
logit_period3_bc <- alpaca::biasCorr(logit_period3)

logit_period3_bc_vcov <- vcov(logit_period3_bc, type = "sandwich", cluster = ~dyad)
logit_period3_bc_vcov_rse <- sqrt(diag(logit_period3_bc_vcov))

summary_logit_period3_bc_rcm <- summary(logit_period3_bc, type = "sandwich", cluster = ~dyad)$cm
logit_period3_bc_rpvalues <- summary_logit_period3_bc_rcm[,4]

logit_period3_bc_apes <- getAPEs(logit_period3_bc)

logit_period3_bc_apes_vcov <- vcov(logit_period3_bc_apes, type = "sandwich", cluster = ~dyad)
logit_period3_bc_apes_vcov_rse <- sqrt(diag(logit_period3_bc_apes_vcov))

summary_logit_period3_bc_apes_rcm <- summary(logit_period3_bc_apes, type = "sandwich", cluster = ~dyad)$cm
logit_period3_bc_apes_rpvalues <- summary_logit_period3_bc_apes_rcm[,4]

mfx_logit_period3 <- data.frame(Coefficient = summary_logit_period3_bc_apes_rcm[,1],
                              SE = summary_logit_period3_bc_apes_rcm[,2],
                              lower = summary_logit_period3_bc_apes_rcm[,1] - summary_logit_period3_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                              upper = summary_logit_period3_bc_apes_rcm[,1] + summary_logit_period3_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                              Model = "MFX Period Logit 3")
mfx_logit_period3 <- data.frame(Variable_Contrast = row.names(mfx_logit_period3), mfx_logit_period3)
rownames(mfx_logit_period3) = NULL

#### G.1.4. Rare UNCLOS negotiations + dyad FE #####

rare_period1 <- glm(onsmdisp ~ ambilex_factor + 
                      ntodel_lag_log + 
                      unclos_status + 
                      petrexpl + 
                      petractive_lag +
                      amdisp_lag + 
                      atdisp_lag +
                      joint_dem +
                      unclos_all_neg +
                      dyad +  
                      ons_cntr + 
                      ons_cntr_sqd + 
                      ons_cntr_cbd, 
                    data = mbm_phi_wo_F0, 
                    family = binomial(link = "logit"), 
                    method = "brglmFit", 
                    pl = TRUE)

summary(rare_period1)

# Get summary:
coeftest(rare_period1, vcov. = vcovCL(rare_period1, cluster = mbm_phi_wo_F0$year))

# Save robust standard errors:
cov_rare_period1 <- vcovCL(rare_period1, 
                           cluster = mbm_phi_wo_F0$year)

robust_se_rare_period1 <- sqrt(diag(cov_rare_period1))

# Margins:
mfx_rare_period1 <- margins(rare_period1, 
                            variables = c("ambilex_factor", "ntodel_lag_log", "unclos_status", "petrexpl", "petractive_lag",
                                          "amdisp_lag", "atdisp_lag", "joint_dem", "unclos_all_neg",
                                          "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                            vcov = cov_rare_period1)

mfx_rare_period1_df <- data.frame(Variable_Contrast = summary(mfx_rare_period1)[,1],
                                  Coefficient = summary(mfx_rare_period1)[,2],
                                  SE = summary(mfx_rare_period1)[,3],
                                  lower = summary(mfx_rare_period1)[,6],
                                  upper = summary(mfx_rare_period1)[,7],
                                  Model = "MFX Period Rare 1")

rownames(mfx_rare_period1_df) <- NULL

#### G.1.5. Rare UNCLOS in force + dyad FE #####

rare_period2 <- glm(onsmdisp ~ ambilex_factor + 
                      ntodel_lag_log + 
                      petrexpl + 
                      petractive_lag +
                      amdisp_lag + 
                      atdisp_lag +
                      joint_dem +
                      unclos_signed + 
                      unclos_in_force +
                      dyad + 
                      ons_cntr + 
                      ons_cntr_sqd + 
                      ons_cntr_cbd, 
                    data = mbm_phi_wo_F0, 
                    family = binomial(link = "logit"), 
                    method = "brglmFit", 
                    pl = TRUE)

# Get summary:
coeftest(rare_period2, vcov. = vcovCL(rare_period2, cluster = mbm_phi_wo_F0$year))

# Save robust standard errors:
cov_rare_period2 <- vcovCL(rare_period2, 
                           cluster = mbm_phi_wo_F0$year)

robust_se_rare_period2 <- sqrt(diag(cov_rare_period2))

# Margins:
mfx_rare_period2 <- margins(rare_period2, 
                            variables = c("ambilex_factor", "ntodel_lag_log", "petrexpl", "petractive_lag",
                                          "amdisp_lag", "atdisp_lag", "joint_dem", "unclos_signed", "unclos_in_force",
                                          "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                            vcov = cov_rare_period2)

mfx_rare_period2_df <- data.frame(Variable_Contrast = summary(mfx_rare_period2)[,1],
                                  Coefficient = summary(mfx_rare_period2)[,2],
                                  SE = summary(mfx_rare_period2)[,3],
                                  lower = summary(mfx_rare_period2)[,6],
                                  upper = summary(mfx_rare_period2)[,7],
                                  Model = "MFX Period Rare 2")

rownames(mfx_rare_period2_df) <- NULL

#### G.1.6. Rare PostCW + dyad FE #####

rare_period3 <- glm(onsmdisp ~ ambilex_factor + 
                      ntodel_lag_log + 
                      unclos_status + 
                      petrexpl + 
                      petractive_lag +
                      amdisp_lag + 
                      atdisp_lag +
                      joint_dem +
                      post_cw +
                      dyad + 
                      ons_cntr + 
                      ons_cntr_sqd + 
                      ons_cntr_cbd, 
                    data = mbm_phi_wo_F0, 
                    family = binomial(link = "logit"), 
                    method = "brglmFit", 
                    pl = TRUE)

summary(rare_period3)

# Get summary:
coeftest(rare_period3, vcov. = vcovCL(rare_period3, cluster = mbm_phi_wo_F0$year))

# Save robust standard errors:
cov_rare_period3 <- vcovCL(rare_period3, 
                           cluster = mbm_phi_wo_F0$year)

robust_se_rare_period3 <- sqrt(diag(cov_rare_period3))

# Margins:
mfx_rare_period3 <- margins(rare_period3, 
                            variables = c("ambilex_factor", "ntodel_lag_log", "unclos_status", "petrexpl", "petractive_lag",
                                          "amdisp_lag", "atdisp_lag", "joint_dem", "post_cw",
                                          "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                            vcov = cov_rare_period3)

mfx_rare_period3_df <- data.frame(Variable_Contrast = summary(mfx_rare_period3)[,1],
                                  Coefficient = summary(mfx_rare_period3)[,2],
                                  SE = summary(mfx_rare_period3)[,3],
                                  lower = summary(mfx_rare_period3)[,6],
                                  upper = summary(mfx_rare_period3)[,7],
                                  Model = "MFX Period Rare 3")

rownames(mfx_rare_period3_df) <- NULL

#### G.1.7. <<< Table OS S5.1 >>> Regression tables for period models ####

texreg(list(logit_period1_bc, rare_period1, 
            logit_period2_bc, rare_period2,
            logit_period3_bc, rare_period3),
       override.se = list(logit_period1_bc_vcov_rse, 
                          robust_se_rare_period1, 
                          logit_period2_bc_vcov_rse, 
                          robust_se_rare_period2,
                          logit_period3_bc_vcov_rse,
                          robust_se_rare_period3),
       override.pvalues = list(logit_period1_bc_rpvalues, NULL,
                               logit_period2_bc_rpvalues, NULL,
                               logit_period3_bc_rpvalues, NULL),
       custom.coef.map = list("ambilex_factorMedium" = "MEDIUM legal uncertainty (ref: LOW)",
                              "ambilex_factorHigh" = "HIGH legal uncertainty (ref: LOW)",
                              "ntodel_lag_log" = "No. of dyadic boundaries to delimit (lagged & logged)",
                              "unclos_statusJoint UNCLOS" = "Joint UNCLOS (ref: Before UNCLOS)", 
                              "unclos_statusNeither UNCLOS" = "Neither UNCLOS (ref: Before UNCLOS)", 
                              "unclos_statusOne state UNCLOS" = "One state UNCLOS (ref: Before UNCLOS)", 
                              "petrexpl" = "Hydrocarbon exploration",
                              "petractive_lag" = "Hydrocarbon activity",
                              "amdisp_lag" = "Ongoing maritime boundary dispute (lagged)", 
                              "atdisp_lag" = "Ongoing related territorial dispute (lagged)",
                              "joint_demJoint democracy" = "One side DEMOCRATIC (ref: joint AUTOCRACY)", 
                              "joint_demOne side democratic" = "Joint DEMOCRACY (ref: joint AUTOCRACY)",
                              "unclos_all_neg" = "UNCLOS negotiations",
                              "unclos_signed" = "UNCLOS signed",
                              "unclos_in_force" = "UNCLOS in force",
                              "post_cw" = "Post-Cold War"),
       booktabs = TRUE,
       dcolumn = TRUE,
       no.margin = TRUE,
       fontsize = "scriptsize",
       include.deviance = FALSE,
       custom.gof.rows = list("Dyad-fixed effects" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                              "Year-fixed effects" = c("No", "No", "No", "No", "No", "No"),
                              "Cubic polynomial for time" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))

#### G.1.8. Marginal effects for period models #### 

margins_periods <- data.frame(rbind(mfx_logit_period1, mfx_logit_period2, mfx_logit_period3)) |>
  filter(Variable_Contrast != "ons_cntr" & 
           Variable_Contrast != "ons_cntr_sqd" & 
           Variable_Contrast != "ons_cntr_cbd")

margins_periods$Variable_Contrast <- factor(margins_periods$Variable_Contrast, 
                                            levels = c("ambilex_factorMedium", 
                                                       "ambilex_factorHigh", 
                                                       "ntodel_lag_log", 
                                                       "unclos_statusJoint UNCLOS", 
                                                       "unclos_statusNeither UNCLOS", 
                                                       "unclos_statusOne state UNCLOS",
                                                       "petrexpl", 
                                                       "petractive_lag", 
                                                       "amdisp_lag", 
                                                       "atdisp_lag", 
                                                       "joint_demOne side democratic",
                                                       "joint_demJoint democracy",
                                                       "unclos_all_neg",
                                                       "unclos_signed",
                                                       "unclos_in_force",
                                                       "post_cw"),
                                            labels = c("MEDIUM legal uncertainty (ref: LOW)",
                                                       "HIGH legal uncertainty (ref: LOW)",
                                                       "No. of dyadic boundaries to delimit (lagged & logged)",
                                                       "Joint UNCLOS (ref: Before UNCLOS)", 
                                                       "Neither UNCLOS (ref: Before UNCLOS)", 
                                                       "One state UNCLOS (ref: Before UNCLOS)", 
                                                       "Hydrocarbon exploration",
                                                       "Hydrocarbon activity",
                                                       "Ongoing maritime boundary dispute (lagged)", 
                                                       "Ongoing related territorial dispute (lagged)",
                                                       "One side DEMOCRATIC (ref: joint AUTOCRACY)", 
                                                       "Joint DEMOCRACY (ref: joint AUTOCRACY)",
                                                       "UNCLOS I, II, III negotiations",
                                                       "1982 UNCLOS signed",
                                                       "1982 UNCLOS in force",
                                                       "Post-Cold War"))

margins_periods$Model <- factor(margins_periods$Model, 
                                levels = c("MFX Period Logit 3",
                                           "MFX Period Logit 2",
                                           "MFX Period Logit 1"), 
                                labels = c("(5) Alternative periods: Post-Cold War",
                                           "(3) Alternative periods: UNCLOS signed or in force",
                                           "(1) Alternative periods: UNCLOS negotiations"))

mep_alt_periods <- ggplot(margins_periods, aes(colour = Model)) + scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2, linewidth = 1) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1,
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, y = Coefficient, 
                      ymin = Coefficient - SE*interval2, 
                      ymax = Coefficient + SE*interval2, 
                      shape = Model), 
                  lwd = 1/2, position = position_dodge(width = 0.5),
                  fill = "white", size = 1) + 
  scale_shape(solid = TRUE) + scale_x_discrete(limits = rev(levels(margins_periods$Variable_Contrast))) + 
  theme_bw(base_size = 30) + 
  coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 3)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 3)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_alt_periods

#### G.1.8.1. <<< OS FIGURE S5.1 >>> Margins with alternative periods, FE Logit ####

ggsave(mep_alt_periods, 
       filename = "os_figure_S5_1.pdf",
       width = 50, 
       height = 40, units = "cm")

# Rare events: 

margins_periods_rare <- data.frame(rbind(mfx_rare_period1_df, 
                                         mfx_rare_period2_df, 
                                         mfx_rare_period3_df)) |>
  filter(Variable_Contrast != "ons_cntr" & 
           Variable_Contrast != "ons_cntr_sqd" & 
           Variable_Contrast != "ons_cntr_cbd")

margins_periods_rare$Variable_Contrast <- factor(margins_periods_rare$Variable_Contrast, 
                                                 levels = c("ambilex_factorMedium", 
                                                            "ambilex_factorHigh", 
                                                            "ntodel_lag_log", 
                                                            "unclos_statusJoint UNCLOS", 
                                                            "unclos_statusNeither UNCLOS", 
                                                            "unclos_statusOne state UNCLOS",
                                                            "petrexpl", 
                                                            "petractive_lag", 
                                                            "amdisp_lag", 
                                                            "atdisp_lag", 
                                                            "joint_demOne side democratic",
                                                            "joint_demJoint democracy",
                                                            "unclos_all_neg",
                                                            "unclos_signed",
                                                            "unclos_in_force",
                                                            "post_cw"),
                                                 labels = c("MEDIUM legal uncertainty (ref: LOW)",
                                                            "HIGH legal uncertainty (ref: LOW)",
                                                            "No. of dyadic boundaries to delimit (lagged & logged)",
                                                            "Joint UNCLOS (ref: Before UNCLOS)", 
                                                            "Neither UNCLOS (ref: Before UNCLOS)", 
                                                            "One state UNCLOS (ref: Before UNCLOS)", 
                                                            "Hydrocarbon exploration",
                                                            "Hydrocarbon activity",
                                                            "Ongoing maritime boundary dispute (lagged)", 
                                                            "Ongoing related territorial dispute (lagged)",
                                                            "One side DEMOCRATIC (ref: joint AUTOCRACY)", 
                                                            "Joint DEMOCRACY (ref: joint AUTOCRACY)",
                                                            "UNCLOS I, II, III negotiations",
                                                            "1982 UNCLOS signed",
                                                            "1982 UNCLOS in force",
                                                            "Post-Cold War"))

margins_periods_rare$Model <- factor(margins_periods_rare$Model, 
                                     levels = c("MFX Period Rare 3",
                                                "MFX Period Rare 2",
                                                "MFX Period Rare 1"), 
                                     labels = c("(6) Alternative periods: Post-Cold War",
                                                "(4) Alternative periods: UNCLOS signed or in force",
                                                "(2) Alternative periods: UNCLOS negotiations"))

mep_alt_periods_rare <- ggplot(margins_periods_rare, aes(colour = Model)) + scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2, linewidth = 1)
mep_alt_periods_rare <- mep_alt_periods_rare + 
  geom_linerange(aes(x = Variable_Contrast, 
                     ymin = Coefficient - SE*interval1,
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Model), 
                  lwd = 1/2, position = position_dodge(width = 0.5), 
                  fill = "white", size = 1) + 
  scale_shape(solid = TRUE) + scale_x_discrete(limits = rev(levels(margins_periods_rare$Variable_Contrast))) + 
  theme_bw(base_size = 30) +
  coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 3)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 3)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_alt_periods_rare

#### G.1.8.1. <<< OS FIGURE S5.2 >>> Margins with alternative periods, Rare Logit ####
ggsave(mep_alt_periods_rare, 
       filename = "os_figure_S5_2.pdf",
       width = 50, 
       height = 40, units = "cm")

#### G.2. Legal void ####

#### G.2.1. Logit FE with nihilex ####

logit_void <- alpaca::feglm(onsmdisp ~ nihilex + 
                            ambilex_factor + 
                            ons_cntr + 
                            ons_cntr_sqd + 
                            ons_cntr_cbd | dyad | dyad, mbm_phi_wo_F0)
summary(logit_void)
logit_void_bc <- alpaca::biasCorr(logit_void)

logit_void_bc_vcov <- vcov(logit_void_bc, type = "sandwich", cluster = ~dyad)
logit_void_bc_vcov_rse <- sqrt(diag(logit_void_bc_vcov))

summary_logit_void_bc_rcm <- summary(logit_void_bc, type = "sandwich", cluster = ~dyad)$cm
logit_void_bc_rpvalues <- summary_logit_void_bc_rcm[,4]

logit_void_bc_apes <- getAPEs(logit_void_bc)

logit_void_bc_apes_vcov <- vcov(logit_void_bc_apes, type = "sandwich", cluster = ~dyad)
logit_void_bc_apes_vcov_rse <- sqrt(diag(logit_void_bc_apes_vcov))

summary_logit_void_bc_apes_rcm <- summary(logit_void_bc_apes, type = "sandwich", cluster = ~dyad)$cm
logit_void_bc_apes_rpvalues <- summary_logit_void_bc_apes_rcm[,4]

mfx_logit_void <- data.frame(Coefficient = summary_logit_void_bc_apes_rcm[,1],
                           SE = summary_logit_void_bc_apes_rcm[,2],
                           lower = summary_logit_void_bc_apes_rcm[,1] - summary_logit_void_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                           upper = summary_logit_void_bc_apes_rcm[,1] + summary_logit_void_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                           Model = "MFX Void")
mfx_logit_void <- data.frame(Variable_Contrast = row.names(mfx_logit_void), mfx_logit_void)
rownames(mfx_logit_void) = NULL

#### G.2.2. Penalized ML with nihilex ####

rare_void <- glm(onsmdisp ~ nihilex + 
                   ambilex_factor + 
                   factor(dyad) + 
                   ons_cntr + 
                   ons_cntr_sqd + 
                   ons_cntr_cbd, 
                 data = mbm_phi_wo_F0, 
                 family = binomial(link = "logit"), 
                 method = "brglmFit", 
                 pl = TRUE)

# Get summary:
coeftest(rare_void, vcov. = vcovCL(rare_void, cluster = mbm_phi_wo_F0$year))

# Save robust standard errors:
cov_rare_void <- vcovCL(rare_void, 
                        cluster = mbm_phi_wo_F0$year)

robust_se_rare_void <- sqrt(diag(cov_rare_void))

# Margins:
mfx_rare_void <- margins(rare_void, 
                         variables = c("ambilex_factor", "nihilex", 
                                       "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                         vcov = cov_rare_void)

mfx_rare_void_df <- data.frame(Variable_Contrast = summary(mfx_rare_void)[,1],
                               Coefficient = summary(mfx_rare_void)[,2],
                               SE = summary(mfx_rare_void)[,3],
                               lower = summary(mfx_rare_void)[,6],
                               upper = summary(mfx_rare_void)[,7],
                               Model = "MFX Void Rare")

rownames(mfx_rare_void_df) <- NULL

# G.2.3. <<< OS TABLE S5.4 >>> Regression table with nihilex ####

texreg(list(logit_void, rare_void),
       custom.coef.map = list("nihilex" = "Legal void",
                              "ambilex_factorMedium" = "MEDIUM legal uncertainty (ref: LOW)",
                              "ambilex_factorHigh" = "HIGH legal uncertainty (ref: LOW)"),
       override.se = list(logit_void_bc_vcov_rse, robust_se_rare_void),
       override.pvalues = list(logit_void_bc_rpvalues, NULL),
       custom.gof.rows = list("Dyad-fixed effects" = c("Yes", "Yes"),
                              "Year-fixed effects" = c("No", "No"),
                              "Cubic polynomial for time" = c("Yes", "Yes")),
       booktabs = TRUE,
       no.margin = TRUE,
       fontsize = "scriptsize",
       dcolumn = TRUE, 
       stars = c(0.001, 0.01, 0.05, 0.10),
       symbol = "+")

#### G.2.4. Marginal effects with nihilex ####
mfx_void <- data.frame(rbind(mfx_logit_void, mfx_rare_void_df)) |>
  filter(Variable_Contrast != "ons_cntr" & 
           Variable_Contrast != "ons_cntr_sqd" & 
           Variable_Contrast != "ons_cntr_cbd")

mfx_void$Variable_Contrast <- factor(mfx_void$Variable_Contrast, 
                                     levels = c("nihilex", 
                                                "ambilex_factorMedium",
                                                "ambilex_factorHigh"),
                                     labels = c("Legal void",
                                                "MEDIUM legal uncertainty (ref: LOW)",
                                                "HIGH legal uncertainty (ref: LOW)"))

mfx_void$Model <- factor(mfx_void$Model, 
                         levels = c("MFX Void Rare",
                                    "MFX Void"), 
                         labels = c("(2) Penalized ML with legal void" ,
                                    "(1) Logit FE with legal void"))

mep_void <- ggplot(mfx_void, aes(colour = Model)) + scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2, linewidth = 1) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1, 
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, 
                 position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Model),
                  lwd = 1/2, 
                  position = position_dodge(width = 0.5), 
                  fill = "white", size = 1) + 
  scale_shape(solid = TRUE) + 
  scale_x_discrete(limits = rev(levels(mfx_void$Variable_Contrast))) + 
  theme_bw(base_size = 30) + 
  coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 2)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 2)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_void

#### G.2.4.1. <<< OS Figure S5.7 >>> Saving marginal effects with nihilex ####
ggsave(mep_void, 
       filename = "os_figure_S5_7.pdf",
       width = 40, 
       height = 20, units = "cm")

#### G.3. Count of disputes ####

#### G.3.1. Yearly data and descriptives ####

yearly_mbm <- mbm_phi |>
  mutate(joint_unclos = dplyr::if_else(unclos_status == "Joint UNCLOS", 1, 0),
         F0 = dplyr::if_else(status_num == "Fully Delimited & Undisputed", 1, 0)) |>
  group_by(year) |>
  dplyr::summarize(ambilex_factor = head(ambilex_factor, 1),
                   number_dyads = length(dyad),
                   n_onsets = sum(onsmdisp),
                   n_F0d = sum(F0),
                   n_joint_unclos = sum(joint_unclos),
                   unclos_all_neg = mean(unclos_all_neg, na.rm = TRUE),
                   n_islands = sum(island_os, na.rm = TRUE),
                   n_cosconfig_adjacency = sum(cosconfig_factor == "Adjacent coasts" | 
                                                 cosconfig_factor == "Both adjacent and opposite coasts"), 
                   n_petractive_lag = sum(petractive_lag, na.rm = TRUE)) |>
  mutate(prop_joint_unclos = n_joint_unclos/number_dyads) |>
  mutate(prop_islands = n_islands/number_dyads) |>
  mutate(prop_adjacency = n_cosconfig_adjacency/number_dyads) |>
  mutate(cum_n_onsets = cumsum(n_onsets)) |>
  ungroup() |>
  mutate(n_onsets_lag = dplyr::lag(n_onsets, n = 1)) |>
  mutate(n_nF0d_lag = dplyr::lag(n_F0d, n = 1),
         delta_n_dyads = number_dyads - dplyr::lag(number_dyads, n = 1))

# Descriptives: 
yearly_mbm <- as.data.frame(yearly_mbm)
stargazer(yearly_mbm,
          keep = c("year", "n_onsets", "ambilex_factor", 
                   "delta_n_dyads",
                   "prop_joint_unclos",
                   "n_onsets_lag",
                   "n_nF0d_lag"),
          no.space = TRUE)

#### G.3.2. Poisson ####

poisson1 <- glm(n_onsets ~ ambilex_factor + 
                  delta_n_dyads + 
                  prop_joint_unclos + 
                  n_onsets_lag + 
                  n_nF0d_lag,
                data = yearly_mbm, 
                family = "poisson") 

summary(poisson1)
poisson1_slopes <- avg_slopes(poisson1)

#### G.3.3. Negative binomial ####

nb1 <- glm.nb(n_onsets ~ ambilex_factor + 
                delta_n_dyads + 
                prop_joint_unclos + 
                n_onsets_lag + 
                n_nF0d_lag,
              data = yearly_mbm)

summary(nb1)
nb1_slopes <- avg_slopes(nb1)

# Likelihood ratio test: 
pchisq(2 * (logLik(nb1) - logLik(poisson1)), df = 1, lower.tail = FALSE)
# Not the case that binomial should clearly be preferred. 

# G.3.4. Regression table and marginal effects ####

# Changing names: 
poisson1_slopes$term  <- c("ambilex_factorHigh", "ambilex_factorMedium", "delta_n_dyads", "n_nF0d_lag", "n_onsets_lag", "prop_joint_unclos")
nb1_slopes$term <- c("ambilex_factorHigh", "ambilex_factorMedium", "delta_n_dyads", "n_nF0d_lag", "n_onsets_lag", "prop_joint_unclos")

#### G.3.4.1. <<< OS TABLE S5.6 >>> Regression output with marginal effects ####
texreg(list(poisson1, poisson1_slopes, nb1, nb1_slopes),
       single.row = FALSE,
       no.margin = TRUE,
       fontsize = "scriptsize",
       include.deviance = TRUE,
       booktabs = TRUE,
       dcolumn = TRUE,
       custom.header = list("Poisson" = 1:2, "Negative Binomial" = 3:4),
       custom.model.names = c("Coefficients", "Marginal effects",
                              "Coefficients", "Marginal effects"),
       custom.coef.map = list("ambilex_factorMedium" = "MEDIUM legal uncertainty (ref: LOW)",
                              "ambilex_factorHigh" = "HIGH legal uncertainty (ref: LOW)",
                              "delta_n_dyads" = "Change in the number of dyads",
                              "prop_joint_unclos" = "Proportion of dyads with joint UNCLOS",
                              "n_onsets_lag" = "Number of onsets in the preceding year",
                              "n_nF0d_lag" = "Cumulative number of fully delimited boundaries",
                              "(Intercept)" = "Constant"))

#### G.4. Legal origin ####

#### G.4.1. Year FE Logit Robustness ####

mbm_phi_wo_F0 <- mbm_phi_wo_F0 |>
  mutate(legal_system_pair = 
           case_when(
             (legal_system.x == "civil" & legal_system.y == "civil") ~ "Both Civil", 
             (legal_system.x == "common" & legal_system.y == "common") ~ "Both Common", 
             (legal_system.x == "islamic" & legal_system.y == "islamic") ~ "Both Islamic", 
             (legal_system.x == "mixed" & legal_system.y == "mixed") ~ "Both Mixed", 
             (legal_system.x == "civil" & legal_system.y == "common") | (legal_system.x == "common" & legal_system.y == "civil") ~ "Civil-Common",
             (legal_system.x == "civil" & legal_system.y == "islamic") | (legal_system.x == "islamic" & legal_system.y == "civil") ~ "Civil-Islamic",
             (legal_system.x == "civil" & legal_system.y == "mixed") | (legal_system.x == "mixed" & legal_system.y == "civil") ~ "Civil-Mixed",
             (legal_system.x == "common" & legal_system.y == "islamic") | (legal_system.x == "islamic" & legal_system.y == "common") ~ "Common-Islamic",
             (legal_system.x == "common" & legal_system.y == "mixed") | (legal_system.x == "mixed" & legal_system.y == "common") ~ "Common-Mixed",
             (legal_system.x == "islamic" & legal_system.y == "mixed") | (legal_system.x == "mixed" & legal_system.y == "islamic") ~ "Islamic-Mixed"),
         legal_system_match = if_else(legal_system.x == legal_system.y, 1, 0)) 

mbm_phi_wo_F0 <- mbm_phi_wo_F0 |>
  mutate(new_independence = case_when(indy_year.x < 1946 & indy_year.y < 1946 ~ "Both before 1946",
                                      (indy_year.x >= 1946 & indy_year.y < 1946) | (indy_year.y >= 1946 & indy_year.x < 1946) ~ "One after 1946", 
                                      indy_year.x >= 1946 & indy_year.y >= 1946 ~ "Both after 1946")) |>
  mutate(new_independence = factor(new_independence, levels = c("Both before 1946",
                                                                "One after 1946",
                                                                "Both after 1946")))

year_fes2_robust1 <- alpaca::feglm(onsmdisp ~ island_os + 
                                     adjacency + 
                                     cseez + 
                                     ntodel_lag_log + 
                                     petrexpl + 
                                     petractive_lag +
                                     catch_tot_lag_logged + 
                                     amdisp_lag + 
                                     atdisp_lag +
                                     joint_dem +
                                     cinc_ratio +
                                     legal_system_match + 
                                     new_independence +
                                     ons_cntr +
                                     ons_cntr_sqd + 
                                     ons_cntr_cbd | year | year, mbm_phi_wo_F0) 

summary(year_fes2_robust1)

year_fes2_robust1_bc <- alpaca::biasCorr(year_fes2_robust1) # Model
year_fes2_robust1_bc_vcov <- vcov(year_fes2_robust1_bc, type = "sandwich", cluster = ~ year)
year_fes2_robust1_bc_vcov_rse <- sqrt(diag(year_fes2_robust1_bc_vcov)) # Robust standard errors

summary_year_fes2_robust1_bc_rcm <- summary(year_fes2_robust1_bc, type = "sandwich", cluster = ~ year)$cm
year_fes2_robust1_bc_rpvalues <- summary_year_fes2_robust1_bc_rcm[,4] # P-values

year_fes2_robust1_bc_apes <- getAPEs(year_fes2_robust1_bc)
year_fes2_robust1_bc_apes_vcov <- vcov(year_fes2_robust1_bc_apes, type = "sandwich", cluster = ~ year)
year_fes2_robust1_bc_apes_vcov_rse <- sqrt(diag(year_fes2_robust1_bc_apes_vcov))

summary_year_fes2_robust1_bc_apes_rcm <- summary(year_fes2_robust1_bc_apes, type = "sandwich", cluster = ~ year)$cm
year_fes2_robust1_bc_apes_rpvalues <- summary_year_fes2_robust1_bc_apes_rcm[,4]

mfx_year_fes2_robust1 <- data.frame(Coefficient = summary_year_fes2_robust1_bc_apes_rcm[,1],
                                    SE = summary_year_fes2_robust1_bc_apes_rcm[,2],
                                    lower = summary_year_fes2_robust1_bc_apes_rcm[,1] - summary_year_fes2_robust1_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                    upper = summary_year_fes2_robust1_bc_apes_rcm[,1] + summary_year_fes2_robust1_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                    Model = "MFX Year FE 2 Robust")

mfx_year_fes2_robust1 <- data.frame(Variable_Contrast = row.names(mfx_year_fes2_robust1), mfx_year_fes2_robust1)
rownames(mfx_year_fes2_robust1) = NULL

#### G.4.2.Year FE Rare Logit Robustness ####

# Model: 
rare_logit_year_fes2_robust1 <- glm(onsmdisp ~ island_os + 
                                      adjacency + 
                                      cseez + 
                                      ntodel_lag_log + 
                                      petrexpl + 
                                      petractive_lag +
                                      catch_tot_lag_logged + 
                                      amdisp_lag + 
                                      atdisp_lag +
                                      joint_dem +
                                      cinc_ratio +
                                      legal_system_match + 
                                      new_independence +
                                      ons_cntr +
                                      ons_cntr_sqd + 
                                      ons_cntr_cbd +
                                      factor(year),
                                    data = mbm_phi_wo_F0, 
                                    family = binomial(link = "logit"), 
                                    method = "brglmFit", 
                                    pl = TRUE)

# Get summary:
coeftest(rare_logit_year_fes2_robust1, vcov. = vcovCL(rare_logit_year_fes2_robust1, cluster = mbm_phi_wo_F0$year))

# Save robust standard errors:
cov_rare_logit_year_fes2_robust1 <- vcovCL(rare_logit_year_fes2_robust1, 
                                           cluster = mbm_phi_wo_F0$year)

robust_se_rare_logit_year_fes2_robust1 <- sqrt(diag(cov_rare_logit_year_fes2_robust1))

# Margins:
mfx_rare_logit_year_fes2_robust1 <- margins(rare_logit_year_fes2_robust1, 
                                            variables = c("island_os", "adjacency", "cseez",
                                                          "ntodel_lag_log", "petrexpl", "petractive_lag", 
                                                          "amdisp_lag", "atdisp_lag", 
                                                          "joint_dem",
                                                          "catch_tot_lag_logged", 
                                                          "cinc_ratio",
                                                          "legal_system_match", 
                                                          "new_independence"), 
                                            vcov = cov_rare_logit_year_fes2_robust1)

mfx_rare_logit_year_fes2_robust1_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_year_fes2_robust1)[,1],
                                                  Coefficient = summary(mfx_rare_logit_year_fes2_robust1)[,2],
                                                  SE = summary(mfx_rare_logit_year_fes2_robust1)[,3],
                                                  lower = summary(mfx_rare_logit_year_fes2_robust1)[,6],
                                                  upper = summary(mfx_rare_logit_year_fes2_robust1)[,7],
                                                  Model = "MFX Rare Logit Year FE 2 Robust")

rownames(mfx_rare_logit_year_fes2_robust1_df) <- NULL

#### G.4.3. Regression table  ####

# Selectiong GOFs for the rare event logit models:

rare_logit_year_fes2_robust1_sle <- texreg::extract(rare_logit_year_fes2_robust1)
rare_logit_year_fes2_robust1_sle@gof.names <- rare_logit_year_fes2_robust1_sle@gof.names[-(1:8)]
rare_logit_year_fes2_robust1_sle@gof.decimal <- rare_logit_year_fes2_robust1_sle@gof.decimal[-(1:8)]
rare_logit_year_fes2_robust1_sle@gof <- rare_logit_year_fes2_robust1_sle@gof[-(1:8)]

#### G.4.3.1. <<< OS Table >>> Legal origin and recent indepencence robustness checks ####

texreg(list(year_fes2_robust1_bc, 
            rare_logit_year_fes2_robust1_sle),
       override.se = list(year_fes2_robust1_bc_vcov_rse, 
                          robust_se_rare_logit_year_fes2_robust1), 
       override.pvalues = list(year_fes2_robust1_bc_rpvalues, 
                               NULL),
       # single.row = TRUE,
       no.margin = TRUE,
       fontsize = "scriptsize",
       custom.model.names = c("FE Logit", "Penalized ML"),
       custom.coef.map = list("island_os" = "Offshore island",
                              "adjacency" = "Adjacency of coasts",
                              "cseez" = "Beyond territorial sea", 
                              "ntodel_lag_log" = "No. of dyadic boundaries to delimit (lagged & logged)",
                              "petrexpl" = "Hydrocarbon exploration",
                              "petractive_lag" = "Hydrocarbon activity",
                              "catch_tot_lag_logged" = "Total catch (lagged & logged)",
                              "cinc_ratio" = "Capability ratio",
                              "amdisp_lag" = "Ongoing maritime boundary dispute (lagged)", 
                              "atdisp_lag" = "Ongoing related territorial dispute (lagged)",
                              "joint_demJoint democracy" = "One side DEMOCRATIC (ref: joint AUTOCRACY)", 
                              "joint_demOne side democratic" = "Joint DEMOCRACY (ref: joint AUTOCRACY)",
                              "legal_system_match" = "Same domestic legal origin",
                              "new_independenceOne after 1946" = "One state independence after 1946 (ref: Both independent before 1946)",
                              "new_independenceBoth after 1946" = "Both states independence after 1946 (ref: Both independent before 1946)"), 
       gof.names = list("", "", "", "", "", "", ""),
       custom.gof.rows = list("Dyad-fixed effects" = c("No", "No"),
                              "Year-fixed effects" = c("Yes", "Yes"),
                              "Cubic polynomial for time" = c("Yes", "Yes")),
       include.deviance = FALSE,
       booktabs = TRUE,
       dcolumn = TRUE, 
       stars = c(0.001, 0.01, 0.05, 0.10),
       symbol = "+")

#### G.4.4. Marginal effects ####

# Margins FE logit and rare events logits:
margins_logit_rare_logit_year_fe_robust1 <- data.frame(rbind(
  mfx_year_fes2_robust1, mfx_rare_logit_year_fes2_robust1_df)) |>
  filter(Variable_Contrast != "ons_cntr" & 
           Variable_Contrast != "ons_cntr_sqd" & 
           Variable_Contrast != "ons_cntr_cbd")

margins_logit_rare_logit_year_fe_robust1$Variable_Contrast <- factor(margins_logit_rare_logit_year_fe_robust1$Variable_Contrast, 
                                                                   levels = c("island_os",
                                                                              "adjacency",
                                                                              "cseez",
                                                                              "ntodel_lag_log", 
                                                                              "petrexpl", 
                                                                              "petractive_lag", 
                                                                              "amdisp_lag", 
                                                                              "atdisp_lag", 
                                                                              "joint_demOne side democratic",
                                                                              "joint_demJoint democracy",
                                                                              "catch_tot_lag_logged", 
                                                                              "cinc_ratio",
                                                                              "legal_system_match",
                                                                              "new_independenceOne after 1946",
                                                                              "new_independenceBoth after 1946"),
                                                                   labels = c("Dyadic uncert.: Offshore island",
                                                                              "Dyadic uncert.: Adjacency of coasts", 
                                                                              "Dyadic uncert.: Beyond territorial sea", 
                                                                              "No. of dyadic boundaries to delimit (lagged & logged)",
                                                                              "Hydrocarbon exploration",
                                                                              "Hydrocarbon activity",
                                                                              "Ongoing maritime boundary dispute (lagged)", 
                                                                              "Ongoing related territorial dispute (lagged)",
                                                                              "One side DEMOCRATIC (ref: joint AUTOCRACY)", 
                                                                              "Joint DEMOCRACY (ref: joint AUTOCRACY)",
                                                                              "Total catch (lagged & logged)",
                                                                              "Capability ratio",
                                                                              "Same domestic legal origin",
                                                                              "One independence post-46 (ref: Both ind. pre-46)",
                                                                              "Both independent post-46 (ref: Both ind. pre-46)"))

margins_logit_rare_logit_year_fe_robust1$Model <- factor(margins_logit_rare_logit_year_fe_robust1$Model, 
                                                       levels = c("MFX Rare Logit Year FE 2 Robust", 
                                                                  "MFX Year FE 2 Robust"), 
                                                       labels = c("(2) Penalized ML FE with controls",
                                                                  "(1) Logit FE with controls"))

mep_origin <- ggplot(margins_logit_rare_logit_year_fe_robust1, aes(colour = Model)) + 
  scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2, linewidth = 1) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1,
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, 
                 position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Model), 
                  lwd = 1/2, 
                  position = position_dodge(width = 0.5), 
                  fill = "white", size = 1) + 
  scale_shape(solid = TRUE) + 
  scale_x_discrete(limits = rev(levels(margins_logit_rare_logit_year_fe_robust1$Variable_Contrast))) + 
  theme_bw(base_size = 30) + 
  coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 2)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 2)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_origin

#### G.4.4.1. <<< OS FIGURE S5.3 >>> Marginal effects for legal origin and recent indy ####
# Saving: 
ggsave(mep_origin, 
       filename = "os_figure_S5_3.pdf",
       width = 55, 
       height = 35, units = "cm")

#### G.5. Two-way FE/DID #### 

# Creating eight-year intervals (because it is a rare event)
# Interval treated: 3 (1969-1976)
# Variables: 
# - The three dyadic uncertainty factors: island_os, adjacency, cseez
# - ambilex_dyad (if any of the three factors present) 

mbm_phi_wo_F0_did <- mbm_phi_wo_F0 |>
  filter(year >= 1953 & year <= 2016) |>
  mutate(eight_y = case_when(
    year >= 1953 & year <= 1960 ~ 1,
    year >= 1961 & year <= 1968 ~ 2,
    year >= 1969 & year <= 1976 ~ 3,
    year >= 1977 & year <= 1984 ~ 4,
    year >= 1985 & year <= 1992 ~ 5,
    year >= 1993 & year <= 2000 ~ 6,
    year >= 2001 & year <= 2008 ~ 7,
    year >= 2009 & year <= 2016 ~ 8,
  )) |>
  mutate(year_treated = if_else(island_os == 1 | adjacency == 1 | cseez == 1, 1969, 5000)) |>
  mutate(ambilex_dyad = if_else(island_os == 1 | adjacency == 1 | cseez == 1, 1, 0)) |>
  mutate(eight_y_treated = if_else(island_os == 1 | adjacency == 1 | cseez == 1, 3, 5000)) |>
  mutate(eight_y_first_treated = if_else(island_os == 1 | adjacency == 1 | cseez == 1, 3, 0)) |>
  as.data.frame()

#### G.5.1. Any dyadic uncertainty factor ####

sunab_did = feols(onsmdisp ~ sunab(eight_y_treated, eight_y) | ambilex_dyad + eight_y, 
                  mbm_phi_wo_F0_did, 
                  cluster = ~ dyad)

summary(sunab_did)

plot_any_dyadic <- ggiplot(sunab_did, 
                           geom_style = 'errorbar',
                           pt.pch = 19, 
                           ci_width = 1,
                           ref.line.par = list(col = 'red', lwd = 1, lty = 2)) + 
  scale_x_continuous(limits = c(-2.5, 5.5),
                     breaks = seq(-2, 5),
                     labels = c("1953-1960", "1961-1968", "1969-1976", 
                                "1977-1984", "1985-1992", "1993-2000",
                                "2001-2008", "2009-2016")) + 
  scale_y_continuous(limits = c(-0.01, 0.03)) + 
  theme_bw(base_size = 15) +
  labs(x = "Years", 
       y = "ATT (Estimate and 95% Conf. Int.)",
       title = "Dyad-specific uncertainty") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot_any_dyadic
#ggsave(filename = "Dynamic DID - Any dyadic factor.pdf",
#       width = 30,
#       height = 15,
#       units = "cm")

#### G.5.2. Offshore islands ####
sunab_did_island = feols(onsmdisp ~ sunab(eight_y_treated, eight_y) | island_os + eight_y, 
                         mbm_phi_wo_F0_did, 
                         cluster = ~ dyad)

summary(sunab_did_island)

plot_island_os <- ggiplot(sunab_did_island, 
                          geom_style = 'errorbar',
                          pt.pch = 19, 
                          ci_width = 1,
                          ref.line.par = list(col = 'red', lwd = 1, lty = 2)) + 
  scale_x_continuous(limits = c(-2.5, 5.5),
                     breaks = seq(-2, 5),
                     labels = c("1953-1960", "1961-1968", "1969-1976", 
                                "1977-1984", "1985-1992", "1993-2000",
                                "2001-2008", "2009-2016")) + 
  scale_y_continuous(limits = c(-0.01, 0.03)) + 
  theme_bw(base_size = 15) +
  labs(x = "Years", 
       y = "ATT (Estimate and 95% Conf. Int.)",
       title = "Offshore islands") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_island_os
#ggsave(filename = "Dynamic DID - Offshore islands.pdf",
#       width = 30,
#       height = 15,
#       units = "cm")

#### G.5.3. CS/EEZ ####
sunab_did_cseez = feols(onsmdisp ~ sunab(eight_y_treated, eight_y) | cseez + eight_y, 
                        mbm_phi_wo_F0_did, 
                        cluster = ~ dyad)

summary(sunab_did_cseez)

plot_cseez <- ggiplot(sunab_did_cseez, 
                      geom_style = 'errorbar',
                      pt.pch = 19, 
                      ci_width = 1,
                      ref.line.par = list(col = 'red', lwd = 1, lty = 2)) + 
  scale_x_continuous(limits = c(-2.5, 5.5),
                     breaks = seq(-2, 5),
                     labels = c("1953-1960", "1961-1968", "1969-1976", 
                                "1977-1984", "1985-1992", "1993-2000",
                                "2001-2008", "2009-2016")) + 
  scale_y_continuous(limits = c(-0.01, 0.03)) + 
  theme_bw(base_size = 15) +
  labs(x = "Years", 
       y = "ATT (Estimate and 95% Conf. Int.)",
       title = "Continental Shelf/EEZ") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_cseez
#ggsave(filename = "Dynamic DID - CS EEZ.pdf",
#       width = 30,
#       height = 15,
#       units = "cm")

#### G.5.4. Adjacency ####
sunab_did_adjacency = feols(onsmdisp ~ sunab(eight_y_treated, eight_y) | adjacency + eight_y, 
                            mbm_phi_wo_F0_did, 
                            cluster = ~ dyad)

summary(sunab_did_adjacency)

plot_adjacency <- ggiplot(sunab_did_adjacency, 
                          geom_style = 'errorbar',
                          pt.pch = 19, 
                          ci_width = 1,
                          ref.line.par = list(col = 'red', lwd = 1, lty = 2)) + 
  scale_x_continuous(limits = c(-2.5, 5.5),
                     breaks = seq(-2, 5),
                     labels = c("1953-1960", "1961-1968", "1969-1976", 
                                "1977-1984", "1985-1992", "1993-2000",
                                "2001-2008", "2009-2016")) + 
  scale_y_continuous(limits = c(-0.01, 0.03)) + 
  theme_bw(base_size = 15) +
  labs(x = "Years", 
       y = "ATT (Estimate and 95% Conf. Int.)",
       title = "Adjacency") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_adjacency
#ggsave(filename = "Dynamic DID - Adjacency.pdf",
#       width = 30,
#       height = 15,
#       units = "cm")

#### G.5.5. <<< OS FIGURE S5.4 >>> Save figure ####
# Combine them using ggpubr

require(grid)   # for the textGrob() function

dids <- ggarrange(plot_any_dyadic + rremove("ylab") + rremove("xlab"), 
                  plot_adjacency + rremove("ylab") + rremove("xlab"), 
                  plot_cseez + rremove("ylab") + rremove("xlab"), 
                  plot_island_os + rremove("ylab") + rremove("xlab"),
                  ncol = 2, nrow = 2) 

annotate_figure(dids, top = textGrob("Dynamic difference-in-differences for dyad-specific uncertainty", 
                                     vjust = 0.45, 
                                     gp = gpar(cex = 2.4)),
                left = textGrob("ATT (Estimate and 95% CI)", 
                                rot = 90, 
                                vjust = 1, 
                                gp = gpar(cex = 2.4)),
                bottom = textGrob("Years", 
                                  gp = gpar(cex = 2.4)))


ggsave(filename = "os_figure_S5_4.pdf",
       width = 30,
       height = 15,
       units = "cm")

#### G.5.6. Regression output (not reported) ####

# Use texreg: 
summary(sunab_did)
texreg(list(sunab_did, 
            sunab_did_adjacency, 
            sunab_did_cseez, 
            sunab_did_island),
       digits = 3,
       custom.model.names = c("Any dyadic uncertainty", "Adjacency", "CS-EEZ", "Offshore islands"),
       custom.coef.map = list("eight_y::-2" = "1953-1960 (ref: 1961-1968)",
                              "eight_y::0" = "1969-1976 (ref: 1961-1968)",
                              "eight_y::1" = "1977-1984 (ref: 1961-1968)",
                              "eight_y::2" = "1985-1992 (ref: 1961-1968)", 
                              "eight_y::3" = "1993-2000 (ref: 1961-1968)", 
                              "eight_y::4" = "2001-2008 (ref: 1961-1968)", 
                              "eight_y::5" = "2009-2016 (ref: 1961-1968)"),
       booktabs = TRUE,
       dcolumn = TRUE, 
       no.margin = TRUE,
       stars = c(0.001, 0.01, 0.05, 0.10),
       symbol = "+")

#### G.6. Regional patterns ####

#### G.6.1. Descriptive analysis ####

# Northern and Western Europe + Baltic Sea
# Mediterranean and Black Sea + Africa
# Central Pacific and East Asia + Indian Ocean and South East Asia + Persian Gulf
# North America + Middle America and Caribbean + South America
# Multi-Region

# Number of onsets per year:

mbm_phi <- mbm_phi |>
  mutate(region_factor_simple = case_when(region_factor %in% c("North America", "South America", "Middle America and Caribbean") ~ "North America, South America, Middle America and Caribbean",
                                          region_factor %in% c("Central Pacific and East Asia", "Indian Ocean and South East Asia") ~ "Central Pacific and East Asia, Indian Ocean and South East Asia",
                                          region_factor %in% c("Africa", "Mediterranean and Black Sea", "Persian Gulf") ~ "Africa, Mediterranean and Black Sea, Persian Gulf", 
                                          region_factor %in% c("Northern and Western Europe", "Baltic Sea") ~ "Northern and Western Europe, Baltic Sea",
                                          region_factor %in% c("Multi-Region") ~ "Multi-Region")) |>
  mutate(region_factor_simple = factor(region_factor_simple, 
                                       levels = c("North America, South America, Middle America and Caribbean", 
                                                  "Africa, Mediterranean and Black Sea, Persian Gulf", 
                                                  "Central Pacific and East Asia, Indian Ocean and South East Asia", 
                                                  "Northern and Western Europe, Baltic Sea", 
                                                  "Multi-Region")))


mbm_phi_wo_F0 <- mbm_phi_wo_F0 |>
  mutate(region_factor_simple = case_when(region_factor %in% c("North America", "South America", "Middle America and Caribbean") ~ "North America, South America, Middle America and Caribbean",
                                          region_factor %in% c("Central Pacific and East Asia", "Indian Ocean and South East Asia") ~ "Central Pacific and East Asia, Indian Ocean and South East Asia",
                                          region_factor %in% c("Africa", "Mediterranean and Black Sea", "Persian Gulf") ~ "Africa, Mediterranean and Black Sea, Persian Gulf", 
                                          region_factor %in% c("Northern and Western Europe", "Baltic Sea") ~ "Northern and Western Europe, Baltic Sea",
                                          region_factor %in% c("Multi-Region") ~ "Multi-Region")) |>
  mutate(region_factor_simple = factor(region_factor_simple, 
                                       levels = c("North America, South America, Middle America and Caribbean", 
                                                  "Africa, Mediterranean and Black Sea, Persian Gulf", 
                                                  "Central Pacific and East Asia, Indian Ocean and South East Asia", 
                                                  "Northern and Western Europe, Baltic Sea", 
                                                  "Multi-Region")))


mbm_phi_ambilex_onsets_region <- mbm_phi |>
  group_by(year, region_factor_simple) |> 
  add_tally(onsmdisp, name = "number_of_onsets") |>
  group_by(year, number_of_onsets, ambilex, region_factor_simple) |>
  summarise(number_of_dyads = n_distinct(dyad)) |>
  mutate(prop_onsets = number_of_onsets/number_of_dyads) 

# Save the annotations:
annotation <- data.frame(
  x = c((1946 + 1960)/2,
        (1960 + 1969)/2,
        (1969 + 1982)/2,
        (1982 + 1993)/2,
        (1993 + 2016)/2),
  y = rep(4.25, 5),
  label = c("LOW", "MEDIUM", "HIGH", "MEDIUM", "LOW"))

# Plotting onsets over time:
onset_plot_region <- ggplot(mbm_phi_ambilex_onsets_region, aes(x = year)) +
  geom_histogram(aes(y = number_of_onsets), stat = "identity", alpha = 0.75) +
  geom_line(aes(y = number_of_dyads/30), linetype = 1, size = 0.75) + 
  scale_y_continuous(name = "Number of disputes", 
                     breaks = seq(0, 5, 1), limits = c(0, 5),
                     sec.axis = sec_axis(trans = ~.*30,   
                                         name = "Number of dyads",
                                         breaks = seq(0, 150, 30))) + 
  scale_x_continuous(name = "Year of dispute onset", breaks = seq(1946, 2016, 10)) +
  theme_bw(base_size = 18) +
  geom_vline(xintercept = c(1960, 1993), linetype = 2, color = "red") + 
  geom_vline(xintercept = c(1969, 1982), linetype = 1, color = "red") +
  geom_label(data=annotation, aes(x = x, y = y, label = label), size = 4) +
  facet_wrap(~ region_factor_simple, nrow = 5)

#### G.6.1.1. <<< OS FIGURE S5.5 >>> Regional patterns of onset ####
# Save onsets over time:
ggsave(onset_plot_region, 
       filename = "os_figure_S5_5.pdf",
       width = 20, 
       height = 24, units = "cm")

#### G.6.2. Control with region ####

#### G.6.2.1. Region 1 ####

mbm_phi_wo_F0_r1 <- mbm_phi_wo_F0 |>
  filter(region_factor_simple == "North America, South America, Middle America and Caribbean")

length(unique(mbm_phi_wo_F0_r1$dyad)) # 101 dyads

#### G.6.2.1.1. Region 1: Dyad FEs ####
rare_logit_fer1 <- glm(onsmdisp ~ ambilex_factor +
                        ons_cntr +
                        ons_cntr_sqd +
                        ons_cntr_cbd + 
                        factor(dyad),
                      data = mbm_phi_wo_F0_r1, 
                      family = binomial(link = "logit"), 
                      method = "brglmFit")

# Get summary:
coeftest(rare_logit_fer1, vcov. = vcovCL(rare_logit_fer1, cluster = mbm_phi_wo_F0_r1$dyad))

# Save robust standard errors:
cov_rare_logit_fer1 <- vcovCL(rare_logit_fer1, cluster = mbm_phi_wo_F0_r1$dyad)

robust_se_rare_logit_fer1 <- sqrt(diag(cov_rare_logit_fer1))

# Margins:
mfx_rare_logit_fer1 <- margins(rare_logit_fer1, 
                              variables = c("ambilex_factor",
                                            "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                              vcov = cov_rare_logit_fer1)

mfx_rare_logit_fer1_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_fer1)[,1],
                                    Coefficient = summary(mfx_rare_logit_fer1)[,2],
                                    SE = summary(mfx_rare_logit_fer1)[,3],
                                    lower = summary(mfx_rare_logit_fer1)[,6],
                                    upper = summary(mfx_rare_logit_fer1)[,7],
                                    Model = "MFX Region 1 Baseline")

rownames(mfx_rare_logit_fer1_df) <- NULL

#### G.6.2.1.2. Region 1: Year FEs ####

# Model: 
rare_logit_year_fes1 <- glm(onsmdisp ~ island_os + 
                              adjacency + 
                              cseez + 
                              ons_cntr +
                              ons_cntr_sqd +
                              ons_cntr_cbd + 
                              factor(year),
                            data = mbm_phi_wo_F0_r1, 
                            family = binomial(link = "logit"), 
                            method = "brglmFit")

# Get summary:
coeftest(rare_logit_year_fes1, vcov. = vcovCL(rare_logit_year_fes1, cluster = mbm_phi_wo_F0_r1$year))

# Save robust standard errors:
cov_rare_logit_year_fes1 <- vcovCL(rare_logit_year_fes1, 
                                   cluster = mbm_phi_wo_F0_r1$year)

robust_se_rare_logit_year_fes1 <- sqrt(diag(cov_rare_logit_year_fes1))

# Margins:
mfx_rare_logit_year_fes1 <- margins(rare_logit_year_fes1, 
                                    variables = c("island_os", "adjacency", "cseez", 
                                                  "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                                    vcov = cov_rare_logit_year_fes1)

mfx_rare_logit_year_fes1_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_year_fes1)[,1],
                                          Coefficient = summary(mfx_rare_logit_year_fes1)[,2],
                                          SE = summary(mfx_rare_logit_year_fes1)[,3],
                                          lower = summary(mfx_rare_logit_year_fes1)[,6],
                                          upper = summary(mfx_rare_logit_year_fes1)[,7],
                                          Model = "MFX Region 1 Dyadic")

rownames(mfx_rare_logit_year_fes1_df) <- NULL

#### G.6.2.2. Region 2 ####

mbm_phi_wo_F0_r2 <- mbm_phi_wo_F0 |>
  filter(region_factor_simple == "Africa, Mediterranean and Black Sea, Persian Gulf")

length(unique(mbm_phi_wo_F0_r2$dyad)) # 119 dyads

#### G.6.2.2.1. Region 2: Dyad FEs ####
rare_logit_fer2 <- glm(onsmdisp ~ ambilex_factor +
                         ons_cntr +
                         ons_cntr_sqd +
                         ons_cntr_cbd + 
                         factor(dyad),
                       data = mbm_phi_wo_F0_r2, 
                       family = binomial(link = "logit"), 
                       method = "brglmFit")

# Get summary:
coeftest(rare_logit_fer2, vcov. = vcovCL(rare_logit_fer2, cluster = mbm_phi_wo_F0_r2$dyad))

# Save robust standard errors:
cov_rare_logit_fer2 <- vcovCL(rare_logit_fer2, cluster = mbm_phi_wo_F0_r2$dyad)

robust_se_rare_logit_fer2 <- sqrt(diag(cov_rare_logit_fer2))

# Margins:
mfx_rare_logit_fer2 <- margins(rare_logit_fer2, 
                               variables = c("ambilex_factor",
                                             "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                               vcov = cov_rare_logit_fer2)

mfx_rare_logit_fer2_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_fer2)[,1],
                                     Coefficient = summary(mfx_rare_logit_fer2)[,2],
                                     SE = summary(mfx_rare_logit_fer2)[,3],
                                     lower = summary(mfx_rare_logit_fer2)[,6],
                                     upper = summary(mfx_rare_logit_fer2)[,7],
                                     Model = "MFX Region 2 Baseline")

rownames(mfx_rare_logit_fer2_df) <- NULL

#### G.6.2.2.2. Region 2: Year FEs ####

# Model: 
rare_logit_year_fes2 <- glm(onsmdisp ~ island_os + 
                              adjacency + 
                              cseez + 
                              ons_cntr +
                              ons_cntr_sqd +
                              ons_cntr_cbd + 
                              factor(year),
                            data = mbm_phi_wo_F0_r2, 
                            family = binomial(link = "logit"), 
                            method = "brglmFit")

# Get summary:
coeftest(rare_logit_year_fes2, vcov. = vcovCL(rare_logit_year_fes2, cluster = mbm_phi_wo_F0_r2$year))

# Save robust standard errors:
cov_rare_logit_year_fes2 <- vcovCL(rare_logit_year_fes2, 
                                   cluster = mbm_phi_wo_F0_r2$year)

robust_se_rare_logit_year_fes2 <- sqrt(diag(cov_rare_logit_year_fes2))

# Margins:
mfx_rare_logit_year_fes2 <- margins(rare_logit_year_fes2, 
                                    variables = c("island_os", "adjacency", "cseez", 
                                                  "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                                    vcov = cov_rare_logit_year_fes2)

mfx_rare_logit_year_fes2_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_year_fes2)[,1],
                                          Coefficient = summary(mfx_rare_logit_year_fes2)[,2],
                                          SE = summary(mfx_rare_logit_year_fes2)[,3],
                                          lower = summary(mfx_rare_logit_year_fes2)[,6],
                                          upper = summary(mfx_rare_logit_year_fes2)[,7],
                                          Model = "MFX Region 2 Dyadic")

rownames(mfx_rare_logit_year_fes2_df) <- NULL

#### G.6.2.3. Region 3 ####

mbm_phi_wo_F0_r3 <- mbm_phi_wo_F0 |>
  filter(region_factor_simple == "Central Pacific and East Asia, Indian Ocean and South East Asia")

length(unique(mbm_phi_wo_F0_r3$dyad)) # 163 dyads

#### G.6.2.3.1. Region 3: Dyad FEs ####
rare_logit_fer3 <- glm(onsmdisp ~ ambilex_factor +
                         ons_cntr +
                         ons_cntr_sqd +
                         ons_cntr_cbd + 
                         factor(dyad),
                       data = mbm_phi_wo_F0_r3, 
                       family = binomial(link = "logit"), 
                       method = "brglmFit")

# Get summary:
coeftest(rare_logit_fer3, vcov. = vcovCL(rare_logit_fer3, cluster = mbm_phi_wo_F0_r3$dyad))

# Save robust standard errors:
cov_rare_logit_fer3 <- vcovCL(rare_logit_fer3, cluster = mbm_phi_wo_F0_r3$dyad)

robust_se_rare_logit_fer3 <- sqrt(diag(cov_rare_logit_fer3))

# Margins:
mfx_rare_logit_fer3 <- margins(rare_logit_fer3, 
                               variables = c("ambilex_factor",
                                             "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                               vcov = cov_rare_logit_fer3)

mfx_rare_logit_fer3_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_fer3)[,1],
                                     Coefficient = summary(mfx_rare_logit_fer3)[,2],
                                     SE = summary(mfx_rare_logit_fer3)[,3],
                                     lower = summary(mfx_rare_logit_fer3)[,6],
                                     upper = summary(mfx_rare_logit_fer3)[,7],
                                     Model = "MFX Region 3 Baseline")

rownames(mfx_rare_logit_fer3_df) <- NULL

#### G.6.2.3.2. Region 3: Year FEs ####

# Model: 
rare_logit_year_fes3 <- glm(onsmdisp ~ island_os + 
                              adjacency + 
                              cseez + 
                              ons_cntr +
                              ons_cntr_sqd +
                              ons_cntr_cbd + 
                              factor(year),
                            data = mbm_phi_wo_F0_r3, 
                            family = binomial(link = "logit"), 
                            method = "brglmFit")

# Get summary:
coeftest(rare_logit_year_fes3, vcov. = vcovCL(rare_logit_year_fes3, cluster = mbm_phi_wo_F0_r3$year))

# Save robust standard errors:
cov_rare_logit_year_fes3 <- vcovCL(rare_logit_year_fes3, 
                                   cluster = mbm_phi_wo_F0_r3$year)

robust_se_rare_logit_year_fes3 <- sqrt(diag(cov_rare_logit_year_fes3))

# Margins:
mfx_rare_logit_year_fes3 <- margins(rare_logit_year_fes3, 
                                    variables = c("island_os", "adjacency", "cseez", 
                                                  "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                                    vcov = cov_rare_logit_year_fes3)

mfx_rare_logit_year_fes3_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_year_fes3)[,1],
                                          Coefficient = summary(mfx_rare_logit_year_fes3)[,2],
                                          SE = summary(mfx_rare_logit_year_fes3)[,3],
                                          lower = summary(mfx_rare_logit_year_fes3)[,6],
                                          upper = summary(mfx_rare_logit_year_fes3)[,7],
                                          Model = "MFX Region 3 Dyadic")

rownames(mfx_rare_logit_year_fes3_df) <- NULL

#### G.6.2.4. Region 4 ####

mbm_phi_wo_F0_r4 <- mbm_phi_wo_F0 |>
  filter(region_factor_simple == "Northern and Western Europe, Baltic Sea")

length(unique(mbm_phi_wo_F0_r4$dyad)) # 33 dyads

#### G.6.2.4.1. Region 4: Dyad FEs ####
rare_logit_fer4 <- glm(onsmdisp ~ ambilex_factor +
                         ons_cntr +
                         ons_cntr_sqd +
                         ons_cntr_cbd + 
                         factor(dyad),
                       data = mbm_phi_wo_F0_r4, 
                       family = binomial(link = "logit"), 
                       method = "brglmFit")

# Get summary:
coeftest(rare_logit_fer4, vcov. = vcovCL(rare_logit_fer4, cluster = mbm_phi_wo_F0_r4$dyad))

# Save robust standard errors:
cov_rare_logit_fer4 <- vcovCL(rare_logit_fer4, cluster = mbm_phi_wo_F0_r4$dyad)

robust_se_rare_logit_fer4 <- sqrt(diag(cov_rare_logit_fer4))

# Margins:
mfx_rare_logit_fer4 <- margins(rare_logit_fer4, 
                               variables = c("ambilex_factor",
                                             "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                               vcov = cov_rare_logit_fer4)

mfx_rare_logit_fer4_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_fer4)[,1],
                                     Coefficient = summary(mfx_rare_logit_fer4)[,2],
                                     SE = summary(mfx_rare_logit_fer4)[,3],
                                     lower = summary(mfx_rare_logit_fer4)[,6],
                                     upper = summary(mfx_rare_logit_fer4)[,7],
                                     Model = "MFX Region 4 Baseline")

rownames(mfx_rare_logit_fer4_df) <- NULL

#### G.6.2.4.2. Region 4: Year FEs ####

# Model: 
rare_logit_year_fes4 <- glm(onsmdisp ~ island_os + 
                              adjacency + 
                              cseez + 
                              ons_cntr +
                              ons_cntr_sqd +
                              ons_cntr_cbd + 
                              factor(year),
                            data = mbm_phi_wo_F0_r4, 
                            family = binomial(link = "logit"), 
                            method = "brglmFit")

# Get summary:
coeftest(rare_logit_year_fes4, vcov. = vcovCL(rare_logit_year_fes4, cluster = mbm_phi_wo_F0_r4$year))

# Save robust standard errors:
cov_rare_logit_year_fes4 <- vcovCL(rare_logit_year_fes4, 
                                   cluster = mbm_phi_wo_F0_r4$year)

robust_se_rare_logit_year_fes4 <- sqrt(diag(cov_rare_logit_year_fes4))

# Margins:
mfx_rare_logit_year_fes4 <- margins(rare_logit_year_fes4, 
                                    variables = c("island_os", "adjacency", "cseez", 
                                                  "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                                    vcov = cov_rare_logit_year_fes4)

mfx_rare_logit_year_fes4_df <- data.frame(Variable_Contrast = summary(mfx_rare_logit_year_fes4)[,1],
                                          Coefficient = summary(mfx_rare_logit_year_fes4)[,2],
                                          SE = summary(mfx_rare_logit_year_fes4)[,3],
                                          lower = summary(mfx_rare_logit_year_fes4)[,6],
                                          upper = summary(mfx_rare_logit_year_fes4)[,7],
                                          Model = "MFX Region 4 Dyadic")

rownames(mfx_rare_logit_year_fes4_df) <- NULL

#### G.6.2.5. Regression output ####

#### G.6.2.5.1 <<< OS TABLE S5.3 >>> Regression output for regional tests ####

texreg(list(rare_logit_fer1, rare_logit_year_fes1,
            rare_logit_fer2, rare_logit_year_fes2,
            rare_logit_fer3, rare_logit_year_fes3,
            rare_logit_fer4, rare_logit_year_fes4),
       override.se = list(robust_se_rare_logit_fer1, robust_se_rare_logit_year_fes1, 
                          robust_se_rare_logit_fer2, robust_se_rare_logit_year_fes2, 
                          robust_se_rare_logit_fer3, robust_se_rare_logit_year_fes3, 
                          robust_se_rare_logit_fer4, robust_se_rare_logit_year_fes4), 
       no.margin = TRUE,
       custom.header = list("Region 1" = 1:2, "Region 2" = 3:4, "Region 3" = 5:6, "Region 4" = 7:8),
       custom.coef.map = list("ambilex_factorMedium" = "MEDIUM legal uncertainty (ref: LOW)",
                              "ambilex_factorHigh" = "HIGH legal uncertainty (ref: LOW)",
                              "island_os" = "Offshore islands",
                              "adjacency" = "Adjacency",        
                              "cseez" = "CS/EEZ"),
       custom.gof.rows = list("Dyad-fixed effects" = c("Yes", "No", "Yes", "No", "Yes", "No", "Yes", "No"),
                              "Year-fixed effects" = c("No", "Yes", "No", "Yes", "No", "Yes", "No", "Yes"),
                              "Cubic polynomial for time" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
       include.deviance = FALSE,
       booktabs = TRUE,
       dcolumn = TRUE, 
       stars = c(0.001, 0.01, 0.05, 0.10),
       symbol = "+")
       
#### G.6.2.6. Marginal effects ####

margins_regions_baseline <- data.frame(rbind(mfx_rare_logit_fer1_df, mfx_rare_logit_fer2_df,
                                             mfx_rare_logit_fer3_df, mfx_rare_logit_fer4_df)) |>
  filter(Variable_Contrast != "ons_cntr" & 
           Variable_Contrast != "ons_cntr_sqd" & 
           Variable_Contrast != "ons_cntr_cbd")

interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

margins_regions_baseline$Variable_Contrast <- factor(margins_regions_baseline$Variable_Contrast, 
                                                           levels = c("ambilex_factorMedium", 
                                                                      "ambilex_factorHigh"),
                                                           labels = c("MEDIUM legal uncertainty (ref: LOW)",
                                                                      "HIGH legal uncertainty (ref: LOW)"))

margins_regions_baseline$Region <- factor(margins_regions_baseline$Model, 
                                               levels = c("MFX Region 4 Baseline", 
                                                          "MFX Region 3 Baseline", 
                                                          "MFX Region 2 Baseline",
                                                          "MFX Region 1 Baseline"), 
                                               labels = c("(4) Northern and Western Europe, Baltic Sea",
                                                          "(3) Central Pacific and East Asia, Indian Ocean and South East Asia" , 
                                                          "(2) Africa, Mediterranean and Black Sea, Persian Gulf",
                                                          "(1) North America, South America, Middle America and Caribbean"))

mep_regions_baseline <- ggplot(margins_regions_baseline, aes(colour = Region)) + 
  scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1,
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast, 
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, 
                      y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Region), 
                  lwd = 1/2, 
                  position = position_dodge(width = 0.5), 
                  fill = "white", size = 1) + 
  scale_shape(solid = TRUE) + scale_x_discrete(limits = rev(levels(margins_regions_baseline$Variable_Contrast))) + 
  theme_bw(base_size = 21) + 
  coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("AME (for models with baseline uncertainty)") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 4)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 4)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_regions_baseline


# Margins logit_fe and rare events logits:
margins_regions_dyadic <- data.frame(rbind(mfx_rare_logit_year_fes1_df, mfx_rare_logit_year_fes2_df,
                                           mfx_rare_logit_year_fes3_df, mfx_rare_logit_year_fes4_df)) |>
  filter(Variable_Contrast != "ons_cntr" & 
           Variable_Contrast != "ons_cntr_sqd" & 
           Variable_Contrast != "ons_cntr_cbd")

margins_regions_dyadic$Variable_Contrast <- factor(margins_regions_dyadic$Variable_Contrast, 
                                                     levels = c("island_os", 
                                                                "adjacency",
                                                                "cseez"),
                                                     labels = c("Offshore islands",
                                                                "Adjacency",
                                                                "Beyond territorial sea"))

margins_regions_dyadic$Region <- factor(margins_regions_dyadic$Model, 
                                         levels = c("MFX Region 4 Dyadic", 
                                                    "MFX Region 3 Dyadic", 
                                                    "MFX Region 2 Dyadic",
                                                    "MFX Region 1 Dyadic"), 
                                         labels = c("(4) Northern and Western Europe, Baltic Sea",
                                                    "(3) Central Pacific and East Asia, Indian Ocean and South East Asia" , 
                                                    "(2) Africa, Mediterranean and Black Sea, Persian Gulf",
                                                    "(1) North America, South America, Middle America and Caribbean"))

mep_regions_dyadic <- ggplot(margins_regions_dyadic, aes(colour = Region)) + 
  scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1, 
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, 
                 position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, 
                      y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Region), 
                  lwd = 1/2, position = position_dodge(width = 0.5), 
                  fill = "white", size = 1) + 
  scale_shape(solid = TRUE) + 
  scale_x_discrete(limits = rev(levels(margins_regions_dyadic$Variable_Contrast))) + 
  theme_bw(base_size = 21) + 
  coord_flip() + 
  ggtitle("") + 
  xlab("") + 
  ylab("AME (for models with dyad-specific uncertainty)") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 4)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 4)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_regions_dyadic

#### G.6.2.6.1 <<< OS FIGURE S5.6 >>> Marginal effects for uncertainty in regions ####

ggarrange(mep_regions_baseline, 
          mep_regions_dyadic,nrow = 2, 
          common.legend = TRUE,
          legend = "bottom")

ggsave(filename = "os_figure_S5_6.pdf",
       width = 30, 
       height = 20, units = "cm")

#### H. PROBING THE MECHANISM ####

#### H.1. Creating the variables ####

mbm_phi_wo_F0 <- mbm_phi_wo_F0 |>
  replace_na(list(spread_policy.x = "No policy", spread_policy.y = "No policy")) |> 
  mutate(TS_match = case_when(TS.x == "No claim" & TS.y == "No claim" ~ "Both without claim",
                              (TS.x != "No claim" & TS.y == "No claim") | (TS.x == "No claim" & TS.y != "No claim") ~ "One without claim",
                              TS.x == "Unspecified" | TS.y == "Unspecified" ~ "At least one Unspecified",
                              TS.x == TS.y ~ "Match", 
                              TS.x %in% c("COORD", "DLM") | TS.y %in% c("COORD", "DLM") ~ "Match", # Assuming these to only be possible when they match
                              TRUE ~ "Mismatch")) |>
  mutate(FZEEZ_match = case_when(FZEEZ.x == "No claim" & FZEEZ.y == "No claim" ~ "Both without claim",
                                 FZEEZ.x == "No claim" | FZEEZ.y == "No claim" ~ "One without claim",
                                 FZEEZ.x %in% c("COORD", "DLM") | FZEEZ.y %in% c("COORD", "DLM") ~ "At least one COORD or DLM",
                                 FZEEZ.x == "Unspecified" | FZEEZ.y == "Unspecified" ~ "At least one Unspecified",
                                 FZEEZ.x == FZEEZ.y ~ "Match",
                                 TRUE ~ "Mismatch")) |>
  mutate(CS_match = case_when(CS.x == "No claim" & CS.y == "No claim" ~ "Both without claim",
                              CS.x == "No claim" | CS.y == "No claim" ~ "One without claim",
                              CS.x %in% c("COORD", "DLM") | CS.y %in% c("COORD", "DLM") ~ "At least one COORD or DLM",
                              CS.x == "Unspecified" | CS.y == "Unspecified" ~ "At least one Unspecified",
                              CS.x == CS.y ~ "Match",
                              CS.x == "CM/200+" & CS.y == "CM/200" | 
                                CS.x == "CM/200" & CS.y == "CM/200+" ~ "Match", 
                              TRUE ~ "Mismatch")) |>
  mutate(Policy_match = case_when(spread_policy.x == "No policy" & spread_policy.y == "No policy" ~ "Both without policy",
                                  (spread_policy.x != "No policy"  & spread_policy.y == "No policy") | (spread_policy.x == 0  & spread_policy.y != 0) ~ "One without policy",
                                  spread_policy.x == spread_policy.y ~ "Match",
                                  spread_policy.x != spread_policy.y ~ "Mismatch")) |>
  mutate(TS_yesmatch = if_else(TS_match == "Match", 1, 0),
         CS_yesmatch = if_else(CS_match == "Match", 1, 0),
         FZEEZ_yesmatch = if_else(FZEEZ_match == "Match", 1, 0),
         Policy_yesmatch = if_else(Policy_match == "Match", 1, 0)) |>
  mutate(all_limits_match = if_else(TS_match == "Match" # 
                                    & FZEEZ_match == "Match", 1, 0)) |>
  mutate(match = case_when(Policy_match == "Match" & all_limits_match == 1 ~ "Both limits and methods match",
                           Policy_match != "Match" & all_limits_match == 1 ~ "Only methods match",
                           Policy_match == "Match" & all_limits_match == 0 ~ "Only limits match",
                           Policy_match != "Match" & all_limits_match == 0 ~ "Neither limits nor methods match")) |>
  arrange(dyad, year) |>
  group_by(dyad) |>
  # to focus on when there are new policies 
  mutate(match_lag = dplyr::lag(match, n = 1)) |>
  mutate(limits_match_methods_match = if_else(match == "Both limits and methods match", 1, 0),
         limits_nomatch_methods_match = if_else(match == "Only methods match", 1, 0),
         limits_match_methods_nomatch = if_else(match == "Only limits match", 1, 0),
         limits_nomatch_methods_nomatch = if_else(match == "Neither limits nor methods match", 1, 0)) 

#### H.2. Descriptive graph of matches over time #### 

mbm_phi_wo_F0 |> group_by(year) |>
  summarize(n_dyads = n_distinct(dyad),
            p_TS_match = sum(TS_match == "Match")/n_dyads,
            p_CS_match = sum(CS_match == "Match")/n_dyads,
            p_FZEEZ_match = sum(FZEEZ_match == "Match")/n_dyads,
            p_Policy_match = sum(Policy_match == "Match")/n_dyads) |>
  pivot_longer(cols = starts_with("p_"),
               names_to = "zone",
               values_to = "p_matches") |>
  mutate(zone = factor(zone, 
                       levels = c("p_TS_match", "p_FZEEZ_match",
                                  "p_CS_match", "p_Policy_match"),
                       labels = c("Territorial sea outer limits", "Fishing zone/EEZ outer limits", 
                                  "Continental shelf (CS) definition", " CS delimitation method"))) |>
  ggplot(aes(x = year)) + 
  geom_line(aes(y = p_matches), stat = "identity", size = 1) + 
  geom_vline(xintercept = c(1969, 1982),
             colour = c("red"),
             linetype = 2,
             size = 1) +
  scale_y_continuous(name = "Proportion of dyad with matching limits or delim. preferences",
                     limits = c(0, 1)) + 
  scale_x_continuous(name = "Year",
                     breaks = seq(1945, 2015, 10)) +
  facet_wrap(~ zone) + 
  theme_bw(base_size = 18) 

#### H.2.1. <<< OS FIGURE S6.1 >>> Save graph for OS ####

ggsave(filename = "os_figure_S6_1.pdf",
       width = 30,
       height = 18,
       units = "cm")

#### H.3. Uncertainty on the probability of matching ####

#### H.3.1. FZ/EEZ ####

logit_FZEEZ_match <- alpaca::feglm(FZEEZ_yesmatch ~ ambilex_factor | 
                                     dyad | 
                                     dyad, 
                                   subset(mbm_phi_wo_F0,
                                          new_FZEEZ.x == 1 | new_FZEEZ.y == 1)) # Only when one or the other side has a new act
summary(logit_FZEEZ_match)

logit_FZEEZ_match_bc <- alpaca::biasCorr(logit_FZEEZ_match) # Model
logit_FZEEZ_match_bc_vcov <- vcov(logit_FZEEZ_match_bc, type = "sandwich", cluster = ~ dyad)
logit_FZEEZ_match_bc_vcov_rse <- sqrt(diag(logit_FZEEZ_match_bc_vcov)) # Robust standard errors

summary_logit_FZEEZ_match_bc_rcm <- summary(logit_FZEEZ_match_bc, type = "sandwich", cluster = ~ dyad)$cm
logit_FZEEZ_match_bc_rpvalues <- summary_logit_FZEEZ_match_bc_rcm[,4] # P-values

logit_FZEEZ_match_bc_apes <- getAPEs(logit_FZEEZ_match_bc)
logit_FZEEZ_match_bc_apes_vcov <- vcov(logit_FZEEZ_match_bc_apes, type = "sandwich", cluster = ~ dyad)
logit_FZEEZ_match_bc_apes_vcov_rse <- sqrt(diag(logit_FZEEZ_match_bc_apes_vcov))

summary_logit_FZEEZ_match_bc_apes_rcm <- summary(logit_FZEEZ_match_bc_apes, type = "sandwich", cluster = ~ dyad)$cm
logit_FZEEZ_match_bc_apes_rpvalues <- summary_logit_FZEEZ_match_bc_apes_rcm[,4]

mfx_logit_FZEEZ_match <- data.frame(Coefficient = summary_logit_FZEEZ_match_bc_apes_rcm[,1],
                                    SE = summary_logit_FZEEZ_match_bc_apes_rcm[,2],
                                    lower = summary_logit_FZEEZ_match_bc_apes_rcm[,1] - summary_logit_FZEEZ_match_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                    upper = summary_logit_FZEEZ_match_bc_apes_rcm[,1] + summary_logit_FZEEZ_match_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                    Model = "Match FZEEZ")

mfx_logit_FZEEZ_match <- data.frame(Variable_Contrast = row.names(mfx_logit_FZEEZ_match), mfx_logit_FZEEZ_match)
rownames(mfx_logit_FZEEZ_match) = NULL

#### H.3.2. CS ####

logit_CS_match <- alpaca::feglm(CS_yesmatch ~ ambilex_factor | 
                                  dyad | 
                                  dyad, 
                                subset(mbm_phi_wo_F0, 
                                       new_CS.x == 1 | new_CS.y == 1))
summary(logit_CS_match)

logit_CS_match_bc <- alpaca::biasCorr(logit_CS_match) # Model
logit_CS_match_bc_vcov <- vcov(logit_CS_match_bc, type = "sandwich", cluster = ~ dyad)
logit_CS_match_bc_vcov_rse <- sqrt(diag(logit_CS_match_bc_vcov)) # Robust standard errors

summary_logit_CS_match_bc_rcm <- summary(logit_CS_match_bc, type = "sandwich", cluster = ~ dyad)$cm
logit_CS_match_bc_rpvalues <- summary_logit_CS_match_bc_rcm[,4] # P-values

logit_CS_match_bc_apes <- getAPEs(logit_CS_match_bc)
logit_CS_match_bc_apes_vcov <- vcov(logit_CS_match_bc_apes, type = "sandwich", cluster = ~ dyad)
logit_CS_match_bc_apes_vcov_rse <- sqrt(diag(logit_CS_match_bc_apes_vcov))

summary_logit_CS_match_bc_apes_rcm <- summary(logit_CS_match_bc_apes, type = "sandwich", cluster = ~ dyad)$cm
logit_CS_match_bc_apes_rpvalues <- summary_logit_CS_match_bc_apes_rcm[,4]

mfx_logit_CS_match <- data.frame(Coefficient = summary_logit_CS_match_bc_apes_rcm[,1],
                                 SE = summary_logit_CS_match_bc_apes_rcm[,2],
                                 lower = summary_logit_CS_match_bc_apes_rcm[,1] - summary_logit_CS_match_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                 upper = summary_logit_CS_match_bc_apes_rcm[,1] + summary_logit_CS_match_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                 Model = "Match CS")

mfx_logit_CS_match <- data.frame(Variable_Contrast = row.names(mfx_logit_CS_match), mfx_logit_CS_match)
rownames(mfx_logit_CS_match) = NULL

#### H.3.3. TS ####

logit_TS_match <- alpaca::feglm(TS_yesmatch ~ ambilex_factor | 
                                  dyad | 
                                  dyad, 
                                subset(mbm_phi_wo_F0, 
                                       new_TS.x == 1 | new_TS.y == 1))
summary(logit_TS_match)

logit_TS_match_bc <- alpaca::biasCorr(logit_TS_match) # Model
logit_TS_match_bc_vcov <- vcov(logit_TS_match_bc, type = "sandwich", cluster = ~ dyad)
logit_TS_match_bc_vcov_rse <- sqrt(diag(logit_TS_match_bc_vcov)) # Robust standard errors

summary_logit_TS_match_bc_rcm <- summary(logit_TS_match_bc, type = "sandwich", cluster = ~ dyad)$cm
logit_TS_match_bc_rpvalues <- summary_logit_TS_match_bc_rcm[,4] # P-values

logit_TS_match_bc_apes <- getAPEs(logit_TS_match_bc)
logit_TS_match_bc_apes_vcov <- vcov(logit_TS_match_bc_apes, type = "sandwich", cluster = ~ dyad)
logit_TS_match_bc_apes_vcov_rse <- sqrt(diag(logit_TS_match_bc_apes_vcov))

summary_logit_TS_match_bc_apes_rcm <- summary(logit_TS_match_bc_apes, type = "sandwich", cluster = ~ dyad)$cm
logit_TS_match_bc_apes_rpvalues <- summary_logit_TS_match_bc_apes_rcm[,4]

mfx_logit_TS_match <- data.frame(Coefficient = summary_logit_TS_match_bc_apes_rcm[,1],
                                 SE = summary_logit_TS_match_bc_apes_rcm[,2],
                                 lower = summary_logit_TS_match_bc_apes_rcm[,1] - summary_logit_TS_match_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                 upper = summary_logit_TS_match_bc_apes_rcm[,1] + summary_logit_TS_match_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                 Model = "Match TS")

mfx_logit_TS_match <- data.frame(Variable_Contrast = row.names(mfx_logit_TS_match), mfx_logit_TS_match)
rownames(mfx_logit_TS_match) = NULL

#### H.3.4. Policy ####

logit_Policy_match <- alpaca::feglm(Policy_yesmatch ~ ambilex_factor | 
                                      dyad | 
                                      dyad, 
                                    subset(mbm_phi_wo_F0,
                                           new_policy.x == 1 | new_policy.y == 1))
summary(logit_Policy_match)

logit_Policy_match_bc <- alpaca::biasCorr(logit_Policy_match) # Model
logit_Policy_match_bc_vcov <- vcov(logit_Policy_match_bc, type = "sandwich", cluster = ~ dyad)
logit_Policy_match_bc_vcov_rse <- sqrt(diag(logit_Policy_match_bc_vcov)) # Robust standard errors

summary_logit_Policy_match_bc_rcm <- summary(logit_Policy_match_bc, type = "sandwich", cluster = ~ dyad)$cm
logit_Policy_match_bc_rpvalues <- summary_logit_Policy_match_bc_rcm[,4] # P-values

logit_Policy_match_bc_apes <- getAPEs(logit_Policy_match_bc)
logit_Policy_match_bc_apes_vcov <- vcov(logit_Policy_match_bc_apes, type = "sandwich", cluster = ~ dyad)
logit_Policy_match_bc_apes_vcov_rse <- sqrt(diag(logit_Policy_match_bc_apes_vcov))

summary_logit_Policy_match_bc_apes_rcm <- summary(logit_Policy_match_bc_apes, type = "sandwich", cluster = ~ dyad)$cm
logit_Policy_match_bc_apes_rpvalues <- summary_logit_Policy_match_bc_apes_rcm[,4]

mfx_logit_Policy_match <- data.frame(Coefficient = summary_logit_Policy_match_bc_apes_rcm[,1],
                                     SE = summary_logit_Policy_match_bc_apes_rcm[,2],
                                     lower = summary_logit_Policy_match_bc_apes_rcm[,1] - summary_logit_Policy_match_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                     upper = summary_logit_Policy_match_bc_apes_rcm[,1] + summary_logit_Policy_match_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                     Model = "Match Policy")

mfx_logit_Policy_match <- data.frame(Variable_Contrast = row.names(mfx_logit_Policy_match), mfx_logit_Policy_match)
rownames(mfx_logit_Policy_match) = NULL

#### H.3.5. <<< Main TABLE 4 >>> Regression output  ####

texreg(list(logit_TS_match, 
            logit_FZEEZ_match, 
            logit_CS_match, 
            logit_Policy_match),
       override.se = list(logit_TS_match_bc_vcov_rse, 
                          logit_FZEEZ_match_bc_vcov_rse, 
                          logit_CS_match_bc_vcov_rse, 
                          logit_Policy_match_bc_vcov_rse), 
       override.pvalues = list(logit_TS_match_bc_rpvalues, 
                               logit_FZEEZ_match_bc_rpvalues, 
                               logit_CS_match_bc_rpvalues, 
                               logit_Policy_match_bc_rpvalues),
       no.margin = TRUE,
       custom.model.names = c("TS Match", "FZEEZ Match", "CS Match", "CS DLM Match"),
       custom.coef.map = list("ambilex_factorMedium" = "MEDIUM legal uncertainty (ref: LOW)",
                              "ambilex_factorHigh" = "HIGH legal uncertainty (ref: LOW)"),
       custom.gof.rows = list("Dyad-fixed effects" = c("Yes", "Yes", "Yes", "Yes"),
                              "Year-fixed effects" = c("No", "No", "No", "No"),
                              "Cubic polynomial for time" = c("No", "No", "No", "No")),
       include.deviance = FALSE,
       booktabs = TRUE,
       dcolumn = TRUE, 
       stars = c(0.001, 0.01, 0.05, 0.10),
       symbol = "+")

#### H.3.6. Marginal effects ####

all_match <- data.frame(rbind(mfx_logit_TS_match,
                              mfx_logit_FZEEZ_match,
                              mfx_logit_CS_match,
                              mfx_logit_Policy_match))

interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

all_match$Variable_Contrast <- factor(all_match$Variable_Contrast, 
                                      levels = c("ambilex_factorMedium",
                                                 "ambilex_factorHigh"),
                                      labels = c("MEDIUM legal uncertainty (ref: LOW)",
                                                 "HIGH legal uncertainty (ref: LOW)"))

all_match$Model <- factor(all_match$Model, 
                          levels = c("Match Policy",
                                     "Match CS",
                                     "Match FZEEZ",
                                     "Match TS"), 
                          labels = c("(4) CS delimitation method",
                                     "(3) Continental shelf (CS)",
                                     "(2) Fishing/Exclusive Economic Zone" ,
                                     "(1) Territorial sea"))

mep_ambilex_on_match <- ggplot(all_match, aes(colour = Model)) + scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2, linewidth = 1) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1,
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, 
                      y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Model), 
                  lwd = 1/2, position = position_dodge(width = 0.5), 
                  fill = "white", size = 1) + scale_shape(solid = TRUE) + 
  scale_x_discrete(limits = rev(levels(all_match$Variable_Contrast))) + 
  theme_bw(base_size = 21) + coord_flip() + 
  ggtitle("DV: Matching limits or delimitation methods") + 
  xlab("") + 
  ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 4)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 4)) +
  theme(legend.position = "bottom", legend.box = "vertical")

mep_ambilex_on_match

#### <<< OS FIGURE S6.2 >>> Marginal effects of uncertainty on matching

ggsave("os_figure_S6_2.pdf",
       width = 30,
       height = 15,
       units = "cm")

#### H.4. Matching (and uncertainty) on the probability of onset ####

#### H.4.1. Matching and baseline uncertainty on onset - Logit ####

logitmatch_on_onset <- alpaca::feglm(onsmdisp ~ match_lag + ambilex_factor + ons_cntr + ons_cntr_sqd + ons_cntr_cbd | 
                                       dyad | dyad, mbm_phi_wo_F0)
summary(logitmatch_on_onset)

logitmatch_on_onset_bc <- alpaca::biasCorr(logitmatch_on_onset) # Model
logitmatch_on_onset_bc_vcov <- vcov(logitmatch_on_onset_bc, type = "sandwich", cluster = ~ dyad)
logitmatch_on_onset_bc_vcov_rse <- sqrt(diag(logitmatch_on_onset_bc_vcov)) # Robust standard errors

summary_logitmatch_on_onset_bc_rcm <- summary(logitmatch_on_onset_bc, type = "sandwich", cluster = ~ dyad)$cm
logitmatch_on_onset_bc_rpvalues <- summary_logitmatch_on_onset_bc_rcm[,4] # P-values

logitmatch_on_onset_bc_apes <- getAPEs(logitmatch_on_onset_bc)
logitmatch_on_onset_bc_apes_vcov <- vcov(logitmatch_on_onset_bc_apes, type = "sandwich", cluster = ~ dyad)
logitmatch_on_onset_bc_apes_vcov_rse <- sqrt(diag(logitmatch_on_onset_bc_apes_vcov))

summary_logitmatch_on_onset_bc_apes_rcm <- summary(logitmatch_on_onset_bc_apes, type = "sandwich", cluster = ~ dyad)$cm
logitmatch_on_onset_bc_apes_rpvalues <- summary_logitmatch_on_onset_bc_apes_rcm[,4]

mfx_logitmatch_on_onset <- data.frame(Coefficient = summary_logitmatch_on_onset_bc_apes_rcm[,1],
                                      SE = summary_logitmatch_on_onset_bc_apes_rcm[,2],
                                      lower = summary_logitmatch_on_onset_bc_apes_rcm[,1] - summary_logitmatch_on_onset_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                      upper = summary_logitmatch_on_onset_bc_apes_rcm[,1] + summary_logitmatch_on_onset_bc_apes_rcm[,2]*(-qnorm((1-0.95)/2)), 
                                      Model = "MFX Logit Match FE")

mfx_logitmatch_on_onset <- data.frame(Variable_Contrast = row.names(mfx_logitmatch_on_onset), mfx_logitmatch_on_onset)
rownames(mfx_logitmatch_on_onset) = NULL

#### H.4.2. Matching and baseline uncertainty on onset - Rare Logit ####

# Model: 

rare_logitmatch_fe1 <- glm(onsmdisp ~ match_lag + ambilex_factor +
                        ons_cntr +
                        ons_cntr_sqd +
                        ons_cntr_cbd + 
                        factor(dyad),
                      data = mbm_phi_wo_F0, 
                      family = binomial(link = "logit"), 
                      method = "brglmFit")

# Get summary:
coeftest(rare_logitmatch_fe1, vcov. = vcovCL(rare_logitmatch_fe1, cluster = mbm_phi_wo_F0$dyad))

# Save robust standard errors:
cov_rare_logitmatch_fe1 <- vcovCL(rare_logitmatch_fe1, cluster = mbm_phi_wo_F0$dyad)

robust_se_rare_logitmatch_fe1 <- sqrt(diag(cov_rare_logitmatch_fe1))

# Margins:

mfx_rare_logitmatch_fe1 <- margins(rare_logitmatch_fe1, 
                              variables = c("ambilex_factor", "match_lag",
                                            "ons_cntr", "ons_cntr_sqd", "ons_cntr_cbd"), 
                              vcov = cov_rare_logitmatch_fe1)

mfx_rare_logitmatch_fe1_df <- data.frame(Variable_Contrast = summary(mfx_rare_logitmatch_fe1)[,1],
                                    Coefficient = summary(mfx_rare_logitmatch_fe1)[,2],
                                    SE = summary(mfx_rare_logitmatch_fe1)[,3],
                                    lower = summary(mfx_rare_logitmatch_fe1)[,6],
                                    upper = summary(mfx_rare_logitmatch_fe1)[,7],
                                    Model = "MFX Rare Logit Match FE")

rownames(mfx_rare_logitmatch_fe1_df) <- NULL

#### H.4.5. Regression output ####

# Selecting GOFs for the rare event logit models:
rare_logitmatch_fe1_sle <- texreg::extract(rare_logitmatch_fe1)
rare_logitmatch_fe1_sle@gof.names <- rare_logitmatch_fe1_sle@gof.names[-(1:8)]
rare_logitmatch_fe1_sle@gof.decimal <- rare_logitmatch_fe1_sle@gof.decimal[-(1:8)]
rare_logitmatch_fe1_sle@gof <- rare_logitmatch_fe1_sle@gof[-(1:8)]

#### H.4.5.1. <<< OS TABLE S6.1 >>> Save regression output ####

texreg(list(logitmatch_on_onset_bc, rare_logitmatch_fe1_sle),
       override.se = list(logitmatch_on_onset_bc_vcov_rse, robust_se_rare_logitmatch_fe1),
       override.pvalues = list(logitmatch_on_onset_bc_rpvalues, NULL),
       custom.header = list("DV: Dispute onset" = 1:2),
       custom.model.names = c("FE Logit", "Penalized ML"),
       custom.coef.map = list("ambilex_factorMedium" = "MEDIUM legal uncertainty (ref: LOW)", 
                              "ambilex_factorHigh" = "HIGH legal uncertainty (ref: LOW)",
                              "match_lagNeither limits nor methods match" = "Neither limits nor methods match (ref: Both match)",
                              "match_lagOnly limits match" = "Only limits match (ref: Both match)",
                              "match_lagOnly methods match" = "Only methods match (ref: Both match)"),
       custom.gof.rows = list("Dyad-fixed effects" = c("Yes", "Yes"),
                              "Year-fixed effects" = c("No", "No"),
                              "Cubic polynomial for time" = c("Yes", "Yes")),
       no.margin = TRUE,
       include.deviance = FALSE,
       booktabs = TRUE,
       dcolumn = TRUE, 
       stars = c(0.001, 0.01, 0.05, 0.10),
       symbol = "+")

#### H.4.6. Marginal effects ####

mfx_match_and_ambilex_on_onset <- rbind(mfx_logitmatch_on_onset, mfx_rare_logitmatch_fe1_df) |>
  filter(Variable_Contrast %in% c("match_lagNeither limits nor methods match",
                                  "match_lagOnly limits match",
                                  "match_lagOnly methods match"))

mfx_match_and_ambilex_on_onset$Variable_Contrast <- factor(mfx_match_and_ambilex_on_onset$Variable_Contrast, 
                                                 levels = c("match_lagNeither limits nor methods match",
                                                            "match_lagOnly limits match",
                                                            "match_lagOnly methods match"),
                                                 labels = c("Neither LIMITS nor METHODS match \n(Ref: Both LIMITS and METHODS match)",
                                                            "LIMITS match, METHODS do not match \n(Ref: Both LIMITS and METHODS match)",
                                                            "METHODS match, LIMITS do not match\n(Ref: Both LIMITS and METHODS match)"))
mfx_match_and_ambilex_on_onset$Model <- factor(mfx_match_and_ambilex_on_onset$Model, 
                                                           levels = c("MFX Logit Match FE", "MFX Rare Logit Match FE"),
                                                           labels = c("FE Logit", "Penalized ML"))


mep_match_on_onset <- ggplot(mfx_match_and_ambilex_on_onset, aes(colour = Model)) + 
  scale_colour_manual(values = my_greys) + 
  geom_hline(yintercept = 0, colour = "red", lty = 2, width = 1) + 
  geom_linerange(aes(x = Variable_Contrast,
                     ymin = Coefficient - SE*interval1,
                     ymax = Coefficient + SE*interval1),
                 lwd = 1, 
                 position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = Variable_Contrast,  
                     ymin = Coefficient - SE*interval2,
                     ymax = Coefficient + SE*interval2),
                 lwd = 1/2, 
                 position = position_dodge(width = 0.5)) + 
  geom_pointrange(aes(x = Variable_Contrast, 
                      y = Coefficient, 
                      ymin = Coefficient - SE*interval2,
                      ymax = Coefficient + SE*interval2, 
                      shape = Model), 
                  lwd = 1/2, position = position_dodge(width = 0.5), 
                  fill = "white",
                  size = 1) + 
  scale_shape(solid = TRUE) + 
  scale_x_discrete(limits = rev(levels(mfx_match_and_ambilex_on_onset$Variable_Contrast))) + 
  scale_y_continuous(limits = c(-0.1, 0.1), breaks = seq(-0.1, 0.1, 0.05)) + 
  theme_bw(base_size = 21)

mep_match_on_onset <- mep_match_on_onset + coord_flip() + 
  ggtitle("DV: Dispute onset") + xlab("") + ylab("Average Marginal Effects") + 
  scale_size(guide = 'none') + 
  guides(colour = guide_legend(reverse = TRUE, nrow = 2)) + 
  scale_shape_manual(values=c(21, 22, 23, 24)) + 
  guides(shape = guide_legend(reverse = TRUE, nrow = 2)) +
  theme(legend.position = "bottom", legend.box = "vertical")

#### H.4.6.1. <<< Main FIGURE 5 >>> Save marginal effects ####

ggsave(filename = "main_figure_5.pdf",
       width = 30,
       height = 12,
       units = "cm")